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This Application is a Continuation-In-Part of Application 
No. 10/011,206, "Empirical Mode Decomposition Apparatus, 



L 10 Method and Article of Manufacture for Analyzing Biological g { S~5^ 

13 Signal and Performing Curve Fitting, filed on December 10, & ^ 
A 2001, which itself claims priority under U.S.C. §12 0 to parent 
;{ application Serial No. 09/282,424 filed on March 31, 1999, 
1 which itself claims priority under U.S.C. §120 to parent / t<?l^ 

\S application Serial No. 08/872, 586 filed on June 10, 1997,^^^ 
] which itself claims priority under 35 U.S.C. §119 (e) to U.S. 
Provisional application Serial No. 60/023,411 filed on August 
14, 1996 and Serial No. 60/023,822 filed on August 12, 1996. 
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This application claims priority under 35 U.S.C. §119 (e) to 
20 U.S. Provisional application Serial No. 60/266,422 filed on 
February 14, 2001. All of the above provisional and non- 
provisional patent applications are hereby incorporated by 
reference . 



2 5 ORIGIN OF INVENTION 

The inventor of the invention described herein is an 
employee of the United States Government. Therefore, the 
invention may be manufactured and used by or for the 
Government for governmental purposes without the payment of 
30 any royalties thereon or therefor. 



COLOR DRAWING NOTIFICATION 

The file of this patent contains at least one drawing 
executed in color. Copies of this patent with color 
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drawing (s) will be provided by the Patent and Trademark Office 
upon request and payment of the necessary fee. 

COPYRIGHT NOTIFICATION 

5 Portions of this patent application contain materials 

that are subject to copyright protection. The copyright owner 
has no objection to the facsimile reproduction by anyone of 
the patent document or the patent disclosure, as it appears in 
the Patent and Trademark Office patent file or records, but 
10 otherwise reserves all copyright rights whatsoever. 
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H BACKGROUND OF THE INVENTION 

UJ This invention generally relates to a signal analysis 

=n method. The results of processing several examples of 

"•is biological and acoustical signals are discussed herein to show 

p the particular utility of the invention in that field and to 

{*? further demonstrate the broad applicability of the invention. 

bi Although the present invention finds utility in 

f:! 

?== processing acoustical signals, it is to be understood that any 
20 signal representative of a real world phenomenon such as a 
signal representative of a physical process including 
electrical, mechanical, biological, acoustical, chemical, 
optical, geophysical or other process (es) may be analyzed and 
thereby more fully understood by applying the invention 
25 thereto. The real world signals to which the invention finds 
utility include a wide variety of real world phenomena such as 
the behavior of a stock market, population growth, traffic 
flow, etc.. Furthermore, the term "real world signal" also 
includes "physical signals" representative of physical 
30 processes such as the electrical, mechanical, biological, 
chemical, acoustical, optical, geophysical process (es) 
mentioned above. 

Although the invention is not limited to a particular 
type of signal processing and includes the full range of real 
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world data representative of processes or phenomena or 
combinations thereof, it is most useful when such real world 
signals are nonlinear and nonstationary . 

5 DESCRIPTION OF RELATED ART 

In the parent application, several examples of 
geophysical data signals representative of earthquakes, ocean 
waves, tsunamis, ocean surface elevation and wind were 
processed to show the invention's wide utility to a broad 
L X0 variety of signal types. The techniques disclosed therein and 
□ elaborated upon herein represent major advances in physical 
Cl signal processing. 

4j Previously, analyzing signals, particularly those having 

Iff nonlinear and/or nonstationary properties, was a difficult 
/is problem confronting many industries. These industries have 

Q harnessed various computer implemented methods to process data 

Fl I 

r7 signals measured or otherwise taken from various processes 
such as electrical, mechanical, optical, biological, and 
f[j chemical processes. Unfortunately, previous methods have not 
20 yielded results which are physically meaningful. 

Among the difficulties found in conventional systems is 
that representing physical processes with physical signals may 
present one or more of the following problems: 
(a) The total data span is too short; 
25 (b) The data are nonstationary; and 

(c) The data represent nonlinear processes. 
Although problems (a) -(c) are separate issues, the first 
two problems are related because a data section shorter than 
the longest time scale of a stationary process can appear to 
30 be nonstationary. Because many physical events are transient, 
the data representative of those events are nonstationary. 
For example, a transient event such as an earthquake will 
produce nonstationary data when measured. Nevertheless, the 
nonstationary character of such data is ignored or the effects 
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assumed to be negligible. This assumption may lead to 
inaccurate results and incorrect interpretation of the 
underlying physics as explained below. 

A variety of techniques have been applied to nonlinear, 
5 nonstationary physical signals. For example, many computer 
implemented methods apply Fourier spectral analysis to examine 
the energy- frequency distribution of such signals. 

Although the Fourier transform that is applied by these 
computer implemented methods is valid under extremely general 
10 conditions, there are some crucial restrictions: the system 

Lsi must be linear, and the data must be strictly periodic or 

# = i 

~: stationary. If these conditions are not met, then the 

'~'-4 resulting spectrum will not make sense physically. 

Ui 

in A common technique for meeting the linearity condition is 

*5j5 to approximate the physical phenomena with at least one linear 

%! 

s system. Although linear approximation is an adequate solution 

ri 

ll\ for some applications, many physical phenomena are highly 

H nonlinear and do not admit a reasonably accurate linear 

[j l 

T-l approximation . 

Furthermore, imperfect probes / sensors and numerical 
schemes may contaminate data representative of the phenomenon. 
For example, the interactions of imperfect probes with a 
perfect linear system can make the final data nonlinear. 

Many recorded physical signals are of finite duration, 
25 nonstationary, and nonlinear because they are derived from 
physical processes that are nonlinear either intrinsically or 
through interactions with imperfect probes or numerical 
schemes. Under these conditions, computer implemented methods 
which apply Fourier spectral analysis are of limited use. For 
30 lack of alternatives, however, such methods still apply 
Fourier spectral analysis to process such data. 

In summary, the indiscriminate use of Fourier spectral 
analysis in these methods and the adoption of the stationarity 
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and linearity assumptions may give inaccurate results some of 
which are described below. 

First, the Fourier spectrum defines uniform harmonic 
components globally. Therefore, the Fourier spectrum needs 
many additional harmonic components to simulate nons tationary 
data that are nonuniform globally. As a result, energy is 
spread over a wide frequency range . 

For example, using a delta function to represent the 
flash of light from a lightning bolt will give a phase-locked 
wide white Fourier spectrum. Here, many Fourier components 
are added to simulate the nonstationary nature of the data in 
the time domain, but their existence diverts energy to a much 
wider frequency domain. Constrained by the conservation of 
energy principle, these spurious harmonics and the wide 
frequency spectrum cannot faithfully represent the true energy 
density of the lighting in the frequency and time space. 

More seriously, the Fourier representation also requires 
the existence of negative light intensity so that the 
components can cancel out one another to give the final delta 
function representing the lightning. Thus, the Fourier 

components might make mathematical sense, but they often do 
not make physical sense when applied. 

Although no physical process can be represented exactly 
by a delta function, some physical data such as the near field 
strong earthquake energy signals are of extremely short 
duration. Such earthquake energy signals almost approach a 
delta function, and they always give artificially wide Fourier 
spectra . 

Second, Fourier spectral analysis uses a linear 
superposition of trigonometric functions to represent the 
data. Therefore, additional harmonic components are required 
to simulate deformed wave profiles. Such deformations, as 
will be shown later, are the direct consequence of nonlinear 
effects. Whenever the form of the data deviates from a pure 



sine or cosine function, the Fourier spectrum will contain 
harmonics. Furthermore, both nons tat ionar i ty and 

nonlinearity can induce spurious harmonic components that 
cause unwanted energy spreading and artificial frequency 
5 smearing in the Fourier spectrum. In other words, the 
nonstationary , stochastic nature of biological data suffers 
from conventional signal processing techniques and makes the 
interpretation of the processed data quite difficult. 

10 Biological Signal Analysis 

lZ According to the above background, there is a need for a 

CJ more accurate signal processing technique that produces 

lj results that are more physically meaningful and readily 

%B understood. Biological signals provide another example of 

in 

•■gjs physical signals in which this invention is applicable. 
J n Parent application serial no. 08/872,586 filed June 10, 1997 
f|| and now issued, U.S. Patent No. 5,983,162, illustrates several 
other types of signals in which this invention is applicable. 
Q Namely, the patent provides specific examples of nonlinear, 

fli 

"20 nonstationary geophysical signals which are very difficult to 
analyze with traditional computer implemented techniques 
including earthquake signals, water wave signals, tsunami 
signals, ocean altitude and ocean circulation signals. 

Many of the aforementioned signal processing problems 

25 exist when biological signals are processed. For example, 
most data in the field of biology are nonstationarily 
stochastic. When conventional tools such as Fourier Analysis 
are applied to such biological data, the result often-times 
obscures the underlying processes. In other words, 

30 conventional Fourier analysis of biological data throws away 
or otherwise obscures valuable information. Thus, the complex 
biological phenomena producing such data cannot be readily 
understood and is, in any event, represented imprecisely. The 
interpretation of the results of such conventionally processed 
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data may, therefore, be quite difficult. The conventional 
techniques also make accurate modelling of the biological 
phenomena very difficult and, sometimes, impossible. 

5 Acoustical Signal Analysis 

The idea of recording and transmitting sounds and 
speeches depends, for example, on the variations of the 
density in the air, and currents in the telephone. The 
crucial element is on the variation, without which the sound 
10 would be a monotonic tone that would not carry any 
Q^ySf-j information. Speech <&&^jg&tt of time varying acoustical 

O signals, which are nonstationary and nonlinear. In fact, for 
h] the acoustical signal to carry any information at all, be it 
yd speech or music notes, there must be a time variation in 

= • = 

-IB amplitude and frequency continuously and, may be, subtly. 

Unfortunately, the tools we have to deal with nonstationary 
fll processes are quite limited; therefore, we are forced to make 

s : 

\ t i all kinds of approximations: As a result, there is a conflict 
y between the human perception and automatic sound processes. 

~s i 

20 Some publications listed, relating to acoustical signal 

analysis, are incorporated by reference: 

Allen, J. B. , 1994: How do humans process and recognize 
speech? IEEE Trans. Speech and Audio Proc . , 2, 567-577. 

25 Banbrook, M, S. McLaughlin, and I. Mann, 1999: Speech 

characterization and synthesis by nonlinear methods. 
IEEE Trans. Speech and Audio Proc, 7, 1-17. 

Billa, J. and A. El-Jaroudi, 1998: An analysis of the effect of 
30 basilar membrance nonlinearities on noise suppression, J. Acoust. 

Soc. Am. , 103, 2691-2705. 

Breen, A. 19 92: Speech synthesis models: A Review. Electron. 
Commun. Eng. J., 19-31, Feb, 1992. 

Carmona, R. , W.L. Hwang, and B. Torresani 1997: 
Characterization of signals by the ridges of their 
wavelet transform, IEEE Trans. Signal Processing, 45 
2586--2S89. 

D' Alessandro, C, V. Darsinos and B. Yegnanarayana , 1998: 
Effectiveness of a periodic and aperiodic decomposition 
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method for analysis of voice sources, IEEE Trans. Speech 
and Audio Proc . , 6, 12-23. 

Furui , S. and M. M. Sondhi, 1991: Advances in Speech 
Signal Processing, Marcel Dekker, New York. 

George, D. E. and E. Salari, 1997: Real-time pitch 
extraction of voiced speech, J. Network and Computer 
Applications, 20, 379-387. 

Hoffmann, R. and C. M. Westendorf, 1997: The development 
of analysis methods for speech recognition. Behavioural 
Processes, 39, 113-125. 

15 Huang, N. E., 2000: A New Method For Nonlinear And 

^ Nonstationary Time Series Analysis: Empirical Mode 

ll\ Decomposition and Hilbert Spectral Analysis. Proceedings 

H SPIE Conference on Wavelet, Orlando, April 2000. 
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£0 Huang , N . E . , Z . Shen , and S . R . Long , M . C . Wu , E . H . 

II Shih, Q. Zheng, C. C. Tung, 

ij and H. H. Liu, 1998: The Empirical Mode Decomposition 

9 - Method and the Hilbert Spectrum for Non-stationary 

Q Time Series Analysis, Proc. Roy. Soc . London, A454, 

f|5 903-995. 

hi Huang, N. E., Z. Shen, R. S. Long, 199 9: A New View of 

(3 Nonlinear Water Waves - The Hilbert Spectrum, Ann. Rev. 

ill Fluid Mech. 31, 417-457. 

30 

Kay, S. M. and R. Sudhaker, 1986: A zero crossing-based 
spectrum analyzer, IEEE Trans. Acoust . Speech, Signal 
Processing, 34, 96- 104. 

35 Kedem, B., 1986: Spectral analysis and discrimination by 

zero-crossing, Proc. IEEE, 74, 1477-1493. 

Kim, D. S., S. Y. Lee and R. M. Kil, 1999: Auditory 
processing of speech signals for robust speech regnition 
40 in real-world noisy environments, IEEE Trans. Speech and 

Audio Proc, 7, 55-69. 

Kumar, A. and S. Mullick, 199 6: Nonlinear dynamical 
analysis of speech, J. Acoust. Soc. Amer . , 100, 615-629, 
45 1996. 

Kumaresan, R., 1998: An inverse signal approach to 
computing the envelope of a real valued signal, IEEE 
Signal Processing Letters, 5, 256-259. 

50 
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Kumaresan, R. and A. Rao, 1999: Model based approach to 
envelope and positive instantaneous frequency estimation 
of signal with speech applications, J. Acoust. Soc . Am., 
105, 1912-1924 . 

Li, T. H, and J. D. Gibson, 1996: Speech analysis and 
segmentation by parametric filtering, IEEE Trans. Speech 
and Audio Proc . , 4, 203-213. 

Loughlin, P. J . , J. W. Pitton, and L . E . Atlas, 1994: 
Construction of positive time- frequency distributions, 
IEEE Trans. Signal Proc., 42, 2697-2705. 

Loughlin, P. J., B. Tacer, 1996: On the amplitude- and frequency- 
modulation decomposition of signals, J. Acoust. Soc. Am., 100, 1594- 
1601 . 

Maragos, P., J. F. Kaiser, and T. F. Quatieri, 1993: 
Energy separation in signal modulations with application 
to speech analysis, IEEE Trans. Signal Processing, 41, 
3024-3051. 

Maragos, P. and A. Potamianos, 1999: Fractal dimensions 
of speech sounds: Computation and application to 
automatic speech recognition, J. Acoust. Soc. Am, 105, 
1925-1932. 

Picinbono, B., 1997: On the instantaneous amplitude and 
phase of signal, IEEE Trans. Signal Processing, 45, 552- 
560 . 

Pinter, Istvan, 1996: Perceptual wavelet-representation 
of speech signals and its application to speech 
enhancement. Computer Speech and language, 10, 1-22. 

Pitton, J. W. , L. E. Atlas, and P. J. Loughlin, 
Applocayionsof positive time- frequency distributions to 
speech processing. IEEE Trans. Speech and Audio Proc., 
2, 554-566. 

Poletti, M. A., 1997: The Homomorphic analytic signal, 
IEEE Trans. Signal Processing, 45, 1943-1953. 

Potamianos, A. and P. Maragos, 1994: A comparison of the 
energy operator and Hilbert transform approach to signal 
and speech demodulation, Signal Processing, 37, 95-120. 

Potamianos, A. and P. Maragos, 1999: Speech analysis and 
synthesis using an AM-FM modulation model, Speech 
Communication, 28, 195-209. 



Quatieri, T. F., T. E. Hanna, and G. C. O'Leary, 1997: 
AM-FM separation using auditory-motivated filters, IEEE 
Trans. Speech and Audio Processing, 5, 465-480. 

Rabiner, L. R. and R. W. Schafer, 1978: Digital 
Processing- of Speech Signals, Prentice hall, Upper 
Saddle River, NJ. 512pp. 

Rabiner, L. R. and B. H. Juang, 1993: Fundamentals of 
Speech Recognition, Prentice Hall, Englewoods, NJ, 507pp. 

Schroeter, J. and M. M . Sondhi, 19 94: Techniques for 
estimating vocal-tract shapes from the speech signal, 
IEEE Trans. Speech and Audio Proc . , 2, 133-150. 

Shower, E. G. and R. Biddulph, 1931: Differnttial pitch 
sensitivity of the ear. J. Acoust. Soc . Am, 3, 275-287. 

Welling, L. and H. Ney, 1998: Formant estimation for 
speech recognition, IEEE Trans. Speech and Audio Proc, 
6, 36-48. 

Yegnanarayana , B., and R. N. J. Veldhuis, 1998: 
Extraction of vocal-tract system characteristics from 
speech signals. IEEE Trans. Speech and Audio Proc., 6, 
313-327 . 

Yegnanarayana, B., C. d'Alessandro and V. Darsinos, 1998: 
An iterative althorithm for decomposition of speech 
signals into periodic and aperiodic components, IEEE 
Trans. Speech and Audio Proc., 6, 1-11. 

Because our emphasis will be on speech analysis, we 

should first examine the principles of Human Speech 

Recognition (HSR) and Automatic (Machine) Speech Recognition 

(ASR) . As summarized in the classical paper by Allen (1994) , 

typical ASR systems start with a front end that transforms the 

speech signal into a feature vector. This processes is mostly 

through spectral analysis over a fixed period of time, within 

such period the speech signal is assumed to be stationary. 

The analysis is strictly on frequency. HSR on the other hand, 

processes information across frequency localized in time. 

Thus, the process is assumed to be nonstationazry. These 

localized speech features are known as the formats. To 

extract features localized in time but across all frequencies 

10 




requires time- frequency analysis. The speech signals are 
divided into a time- frequency continuum of formants by the 
cochlea of the ear, which, by its mechanical properties, is a 
very poor frequency discriminator. Yet, classic experiments 
by Shower and Biddulph (1931) has shown that the ear can 
detect frequency difference as small as 3 Hz near a signal of 
1000 Hz. Such perceptual acuity for pitch seems to violate 
the 'uncertainty relation' for stationary processes, where _f 

t __, with _£ as the standard deviation of the frequency, 

and _t as the given time period. Clearly, using the 
traditional method based on the stationary assumption cannot 
explain the HSR processes. Therefore, for an accurate sound 
perception, we need a different paradigm of sound signal 
analysis, one with time- frequency analysis. This is crucial 
for speech recognition. It is also the cause of the long 
standing difficulty in Speech recognition as stated in the 
classic book by Rabiner and Juang (1993) : ''Although there is 
solid basis for the linguist description of sounds and a good 
understanding of the associates acoustics of sound production 
there is, at best, a tenuous relationship between a given 
linguistic sound and a repeatable, reliable, measurable set o 
acoustic parameters." To overcome this difficulty, we need a 
new signal analysis method. 

To deal with the nonstationary properties of speech, 
various methods have been employed (see, for example, the 
classical book by Rabiner and Schafer (1978) , and more recent 
developments by Furui and Sonfhi (1992) , Hoffman and 
Westendorf (1997)), that include spectral analysis, filter- 
bank, zero-crossing, pattern recognition dynamical 
programming, linear prediction, statistic methods, and neural 
network. In many of these approaches, there lies a tactic 
assumption that the speech can be treated as locally 
stationary. Although great progress has been made, in speech 
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recognition, the locally stationary assumption has rendered 
speech synthesis to bear a flat, wooden and artificial tone. 

In addition to speech recognition, there are a wide 
variety of speech communications application that will depend 
5 on detailed signal analysis, such as digital transmission and 
storage of speech signals, speech synthesis, speaker 
verification and identification, and the enhancement of signal 
quality. The techniques used in all of these applications are 
not strictly limited to speeches alone, the same approach can 
10 also be applied to musical scores and recordings, and even to 
t*J machine condition monitoring. 

□ So far, we have only touched on the problem of 

l 5 | nonstationary properties of the acoustical signal. A more 

important part is the nonlinear properties of the acoustical 

in 

SJ|5 * signal in general and speech in particular, which has mostly 
|U been neglected up until very recently (see, for example, 
fll Kumar, 1996) . As it is well known that the sound in speech is 
i a l generated by three mechanisms: 

P 

20 (1) The vibration of the vocal chords, 

(2) The friction of air through construction of the 

vocal tract, and 
(3) The explosion of a sudden release of the air from 
complete closure of the vocal tract. 

25 

Of the three mechanisms, the vowels are generated by the 
vibration of the vocal cord with unrestricted passage of air. 
Such sound can be generated for indefinite length as long as 
the lung can supply the airflow. Therefore, the vowels are 
3 0 the only sounds that could be approximated locally as 

stationary, but only locally. Unfortunately, in our speech, #j 
vowels are not the information carrying sounds, the oonoonant»gr 
are. The consotant's varying of frequency is a necessity for 
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transmitting information. For example, if we write the 
sentence , 

Do you understand what John just said. 

5 

as , 

D- y- - -nd-rst-nd wh-t J-hn -jst s- -d. 

10 Most people can still figure out the meaning of the 

M utterance. On the other hand, if we would write the sentence 
as 

H 
Lil 

u! -o - ou u- -e- - -a- - - - a- -o- - u- -ai- . 

* No one would be able to decipher what is its meaning at 

rt I all. Consequently, in speech recognition, the precise 

^" analysis of the consonant signals is crucial. 

UJ 
CD 
rn 

s 20 Consonants are highly transient and the generating 

mechanisms are all nonlinear. Such nonlinearity has been 
known for a long time (See, for, example, Rabiner and Schafer, 
1978) , but only recently has there been any investigations 
(see, for example, Billa and El-Jaroudi, 1998; and Maragos and 

25 Potamianos, 1999) . Although the locally stationary assumption 
can be used for vowels better than the consonants, we have to 
examine the consonants, for, as we have seen, consonants 
contain far more information than vowels. This creates 
obvious difficulties for the present approach in speech signal 

30 analysis using methods based on linear and stationary 
assumptions such as filter bank or spectrogram etc. 
Considering the generation mechanisms, we can immediately find 
that mechanisms involved in- generating the consonants are all 
nonlinear. Therefore, we need a signal analysis method that 
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is not only applicable to nonstationary but also to nonlinear 
processes. The Empirical Mode Decomposition (EMD) is the only 
method known for the task. Previously, an application of the 
EMD been the Hilbert-Huang Transform (HHT) , as described in 
5 previous patent applications and publications. In speech 

communication, we process information frequency and amplitude 
locally and instantaneously. The natural way is not to wait 
for a long string of data and then to process it spectrally. 
Built in with an intermittent test option, we can use EMD as a 
10 perfect tool to construct a filter bank, but the filter is in 
p temporal space rather than in frequency space. This temporal 
O space filtering retains the full nonlinearity without 



resorting to the spurious harmonics. It also eliminates the 



distortion that could be introduced through Fourier based 
Hjjs frequency filtering. Therefore, EMD is a natural choice. The 
JL product of EMD is actually the formant components without the 
fU difficulties of frequency resolution in any range. 



Speech signal analysis is the most fundamental 



^ requirement for speech synthesis, speaker verification and 

ru 

20 identification, speech recognition, and enhancement and 

restoration of speech record. The basic technique can also be 
applied to processing music 



SUMMARY OF THE INVENTION 

25 An object of the present invention is to solve the above- 

mentioned problems in conventional signal analysis techniques . 

Another object of the present invention is to provide 
further examples of physical signal processing thereby further 
demonstrating the broad applicability of the invention to a 

30 wide array of physical signals which include acoustical 
signals . 

Another object is to provide a technique of distilling a 
physical signal to the point at which the signal can be 
represented with an analytic function. 
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To achieve these objects, the invention employs a 
computer implemented Empirical Mode Decomposition method which 
decomposes physical signals representative of a physical 
phenomenon into components. These components are designated 
5 as Intrinsic Mode Functions (IMFs) and are indicative of 
intrinsic oscillatory modes in the physical phenomenon. 

Contrary to almost all the previous methods, this new 
computer implemented method is intuitive, direct, a 
posteriori , and adaptive, with the basis of the decomposition 
10 based on and derived from the physical signal. The bases so 
derived have no close analytic expressions, and they can only 
be numerically approximated in a specially programmed computer 
by utilizing the inventive methods disclosed herein. 

More specifically, the general method of the invention 

: : : 

includes two main components or steps to analyze the physical 
"L signal without suffering the problems associated with computer 
rlJ implemented Fourier analysis, namely inaccurate interpretation 
\ d of the underlying physics or biology caused in part by energy 
w spreading and frequency smearing in the Fourier spectrum. 
20 The first step is to process the data with the Empirical 

Mode Decomposition ( EMD) method, with which the data are 
decomposed into a number of Intrinsic Mode Function (IMF) 
components. In this way, the signal will be expanded by 

using a basis that is adaptively derived from the signal 
25 itself. 

The second step of the general method of the present 
invention is to apply the Hilbert Transform to the decomposed 
IMF's and construct an energy-frequency- time distribution, 
designated as the Hilbert Spectrum, from which occurrence of 
30 physical events at corresponding times (time localities) will 
be preserved. There is also no close analytic form for the 
Hilbert Spectrum. As explained below, the invention avoids 
this problem by storing numerical approximations in the 
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specially programmed computer by utilizing the inventive 
method. 

The invention also utilizes instantaneous frequency and 
energy to analyze the physical phenomenon rather than the 
5 global frequency and energy utilized by computer implemented 
Fourier spectral analysis. 

Furthermore, a computer implementing the invention, e.g., 
via executing a program in software, to decompose physical 
signals into intrinsic mode functions with EMD and generate a 
10 Hilbert spectrum is also disclosed. Because of the lack of 
ll % close form analytic expression of either the basis functions 

H and the final Hilbert spectrum; computer implementation of the 

%i 

hi inventive methods is an important part of the overall method, 

■ f\ 

*;j however other implementation of the invention may be done. 

HB Still further, the invention may take the form of an 

n article of manufacture. More specifically, the article of 
manufacture is a computer-usable medium, including a computer- 

: 

[=j readable program code embodied therein wherein the computer- 

z;\ readable code causes a computer to execute the inventive 



20 method. 

Once the IMF * s are generated, the invention can then 



eliminates undesired IMF's and thereby generates a filtered 
25 signal from which it is possible to perform a curve fitting 
process. In this way, it is possible to arrive at an analytic 
function which accurately represents the physically important 
components of the original signal. 



30 spectrum generated from IMF decomposed from a first signal may 
be compared with another Hilbert Spectrum to determine 
identification of the first signal. 

Further scope of applicability of the present invention 
will become apparent from the detailed description given 



produce a distilled or otherwise filtered version of the 
original physical signal. This distillation process 



This invention discloses a method in which Hilbert 
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hereinafter. However, it should be understood that the 
detailed description and specific examples, while indicating 
preferred embodiments of the invention, are given by way of 
illustration only, since various changes and modifications 
5 within the spirit and scope of the invention will become 
apparent to those skilled in the art from this detailed 
description. Furthermore, all the mathematic expressions are 
used as a short hand to express the inventive ideas clearly 
and are not limitative of the claimed invention. 

10 

BRIEF DESCRIPTION OF DRAWINGS 

Q The present invention will become more fully understood 

=1 from the detailed description given hereinbelow and the 

HI 

%il accompanying drawings which are given by way of illustration 

i n 

ie only, and thus are not limitative of the present invention, 

; and wherein: 

flj Figure 1(a) is a high-level flowchart describing the 

r* overall inventive method which may be implemented on the 

□ computer system shown in Figure 2; 

- : : 
I i.| 

: 2 r o Figure 1(b) is a high-level flowchart describing the 

Sifting Process which may be implemented on the computer 
system shown in Figure 2; 

Figure 1(c) is a continuation of the high-level flowchart 
in Figure 1 (b) describing the Sifting Process which may be 
25 implemented on the computer system shown in Figure 2 ; 

Figure 1(d) is a high-level flowchart describing EMD 
signal filtering and curve fitting which may be implemented on 
the computer system shown in Figure 2 ; 

Figure 2 is a high-level block diagram of a computer 
30 system which may be programmed with the inventive with the 
result being a special purpose computer; 

Figure 3(a) shows wind speed data in the form of a graph 
plotting wind speed as a function of time for explaining the 
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computer implemented Empirical Mode Decomposition method of 
the invention; 

Figure 3(b) is a graph illustrating the upper envelope, 
lower envelope, envelope mean and original wind speed data 
5 which are utilized to explain the computer implemented 
Empirical Mode Decomposition method of the invention; 

Figures 3(c) -(e) are graphs of the first, second and 
third component signals hi, hll, hl2 , respectively which are 
generated by the Sifting Process of the invention; 
10 Figure 3(f) is a graph of the first intrinsic mode 

III function component which is generated by the Sifting Process 
w of the invention; 

hj Figure 3 (g) is a graph of data with intermittency for 

?« illustrating an optional intermittency test of the invention; 
"ife Figures 3(h)-(j) are graphs of the first, second, and 

r-\ third intrinsic mode functions when the Sifting Process is 
f*' applied to the data of Figure 3 (g) without applying the 

s 

Uj intermittency test option; 

\l\ Figures 3(k)-(m) are graphs of the first, second, and 

20 third intrinsic mode functions when the Sifting Process is 
applied to the data of Figure 3 (g) which applies the 
intermittency test option. 

Figure 4(a) is a graph of a wind speed signal 
which is for explaining the computer implemented Empirical 
25 Mode Decomposition method of the invention; 

Figures 4 (b) - (k) show the wind speed signal and the nine 
intrinsic mode functions which are extracted from the wind 
speed signal by the EMD method of the invention; 

Figures 5(a)-(j) are a series of graphs illustrating the 
30 successive reconstruction of the original wind speed data from 
the intrinsic mode functions; 

Figure 6(a) is the Hilbert Spectrum generated by the 
invention from the wind speed data of Figure 4(a) ; 
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Figure 6 (b) is the conventional Morlet Wavelet spectrum 
generated f rom the wind speed data of Figure 4(a); 

Figure 6(c) shows the Hilbert Spectrum of Figure 6(a) 
after smoothing by a 15x15 weighted Gaussian smoothing filter; 
5 Figure 7 is a comparison of the marginal Hilbert spectrum 

(solid line) and the Fourier spectrum (dotted line) which were 
generated from the wind speed signal of Figure 4(a); 

Figure 8(a) is a graph illustrating the Degree of 
Stationarity and Degree of Statistical Stationarity which were 
10 generated from the wind speed signal of Figure 4(a) with time 
averages of 10, 50, 100 and 300; Figures 8(b) and (c) are 

13 sections of the wind speed data that was used by the invention 

%J 

y to produce the Degree of Stationarity shown in Figure 8(a); 

. ?% 

yd Figures 9 (a) - (d) , (f) and (g) are graphs of blood 

yi 

Nj5 pressure data taken from the pulmonary artery of an normal, 

2 

PI active rat which provide examples of biological data that may 
rt^ be processed by the invention; 

jb=s 

u'l Figure 9 (e) shows an envelope linking the systolic 

^ pressure extrema values for explaining the concepts of the 
20 invention; 

Figure 10(a) shows a conventional Fourier Spectrum 
(energy versus frequency) of the blood pressure data from 
Figure 9(b) for illustrating advantages of the invention; 

Figures 10(b) -(c) show conventional Fourier Spectrums 
25 (energy versus frequency) of the blood pressure data from 

Figures 9(c) -(d) for further illustrating advantages of the 
invention; 

Figure 10(d) shows a conventional Fourier Spectrum (time 
versus frequency) of the blood pressure data from Figure 9 (b) 
30 for illustrating advantages of the invention; 

Figures 10(e) -(f) show conventional three-dimensional 
Fourier Spectrum (amplitude of spectrum as a function of 
frequency in every 1 -minute window on the time- frequency plane 
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in linear scales) of the blood pressure data from Figure 9(b) 
for illustrating advantages of the invention; 

Figure 10(g) is a combined graph directly comparing a 
conventional Fourier Spectrum and a marginal Hilbert Spectrum 
according to the invention calculated from the data of Figure 
9(c) ; 

Figures 11 (a) -(h) are graphs of the first through eighth 
intrinsic mode functions which are extracted from the blood 
pressure signal of Figure 9(c) by the EMD method of the 
invention; 

Figures 12 (a) -(h) are graphs of the first through eighth 
intrinsic mode functions which are extracted from the blood 
pressure signal of Figure 9(d) by the EMD method of the 
invention; 

Figures 13 (a) and (b) are reconstructions of the blood 
pressure signal of Figure 9(d) based on subsets of the 
intrinsic mode functions; 

Figure 13 (c) is another reconstructions of the blood 
pressure signal of Figure 9(d) based on a different subset of 
the intrinsic mode functions plotted together with the 
original signal (dotted line) of Figures 9(d); 

Figure 14(a) is a Hilbert Spectrum of the Figure 9(d) 
blood pressure signal calculated according to the invention; 

Figure 14(b) is a Hilbert Spectrum of the Figure 9(c) 
blood pressure signal calculated according to the invention; 

Figures 15 (a) -(d) are graphs of pulmonary blood pressure 
signals in response to step changes in oxygen concentration in 
the breathing gas; 

Figures 15 (e) - (h) illustrate the inventive sifting 
process as it is applied to the data of Figures 15 (b) ; 

Figures 16 (a) - (p) are graphs of the first through 
sixteenth intrinsic mode functions which are extracted from 
the blood pressure signal of Figure 15(a) by the EMD method of 
the invention; 
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Figures 17 (a) -(f) are mean trends of pulmonary arterial 
blood pressure which are computed according to the invention; 

Figures 18 (a) -(c) are analytic functions derived by the 
invention and representing the indicial response of pulmonary 
5 arterial blood pressure to a step decrease in oxygen 
concentration from 2 0.9 to 10.0%; 

Figures 19 (a) -(c) are analytic functions derived by the 
invention and representing the indicial response of pulmonary 
arterial blood pressure to a step increase in oxygen 
10 concentration from 10.0 to 20.9%; 

(3 Figures 20 (a) -(d) are oscillations about the mean trend 

M: as defined by the invention for k= 1,2, 4 and 6, respectively; 
UJ Figure 21 shows the Hilbert energy spectrum E k (t) 

i |_ s 

{11 according to the invention which is calculated from the 
'"is instantaneous amplitude spectrum of the oscillations about the 
q mean X k (t) with k=6 (the data from Figure 20(d)); 
? y Figures 22 (a) - (h) are graphs of the first through 

Ui sixteenth intrinsic mode functions which are extracted from 
fl": the Hilbert energy spectrum E k (t) of Figure 21 by the EMD 
20 method of the invention; 

Figures 23 (a) -(b) show the three dimensional (amplitude- 
frequency- time ) and two-dimensional (contour of amplitude on 
the frequency-time plane) plot of the Hilbert Spectrum (HHT) 
taken from the data of Figure 20(d); 
25 Figure 24 shows a conventional two dimensional Fourier 

Spectrum (FFT) of the pressure signal in 1-minute segments 
under the assumption that the process is stationary in each 
segment ; 

Figure 25 is a graph of heart pulse interval versus time 
30 taken from a human with sleep apnea and including both a 
normal and abnormal (apnea) data section; 

Figure 2 6 is a graph of a normal heart pulse interval 
versus time; 
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Figures 27 (a) -(h) are graphs of the first through seventh 
intrinsic mode functions and the residue which are extracted 
from the normal heart pulse interval data of Figure 2 6 by the 
EMD method of the invention; 
5 Figure 28 is a blow up of Figure 26; 

Figure 2 9 is a Hilbert Spectrum of the Figure 2 6 normal 
heart pulse interval data calculated according to the 
invention; 

Figure 3 0 is a graph of an abnormal heart pulse interval 
10 versus time; 

H Figures 31 (a) -(h) are graphs of the first through seventh 

ill 

l:i intrinsic mode functions and the residue which are extracted 
= 1 from the abnormal heart pulse interval data of Figure 2 9 by 
V3 the EMD method of the invention; 

m 

gj5 Figure 3 2 is a blow up of Figure 30; 

L a Figure 3 3 is a Hilbert Spectrum of the Figure 3 0 abnormal 

ill 

til heart pulse interval data calculated according to the 
j ai invention; 

P Figure 34 is heart pulse rate data for 45 minute interval 

f|j 

20 during which the patient suffered an epileptic seizure (at 
time index 0 ) ; 

Figures 35 (a) -(i) are graphs of the first through ninth 
intrinsic mode functions which are extracted from the heart 
pulse interval data of Figure 3 4 by the EMD method of the 
25 invention; 

Figures 36 (a) -(c) are the first intrinsic mode function , 
the original heart pulse interval data, and the third 
intrinsic mode function of the Figure 34 epileptic seizure 
data plotted on an expanded scale, respectively; 
30 Figure 37 is a Hilbert Spectrum of the Figure 34 

epileptic seizure, heart pulse interval data calculated 
according to the invention; and 
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Figure 3 8 is a conventional Wavelet Spectrum of the 
Figure 34 epileptic seizure, heart pulse interval data for 
comparison with the inventive result of Figure 37. 

Figure 3 9 is a block diagram describing the inventive 
5 method of speech analysis. 

Figure 40 is a block diagram describing the inventive 
method of speech synthesis. 

Figure 41 is a block diagram describing the inventive 
method of speaker identification. 
10 Figure 42 is a block diagram describing the inventive 

method of speech recognition. 

Figure 43 is a block diagram describing the inventive 
iy method of sound quality enhancement. 

ul Figure 44 is a block diagram describing the inventive 

Ui 

^a|5 method of machine health monitoring. 

]L Figure 45(a) shows acoustical signal data recorded from a 

fU male speaker, saying 'Halloo.' 

[ a j Figure 45 (b) shows acoustical signal data recorded from a 

female speaker, saying 'Halloo.' 

II 

20 Figure 46 (a) - (b) show the EMD decomposed IMF components 

of Figure 45 (a) - (b) . 

Figure 47 (a) - (b) show the Hilbert Spectra of Figure 
46(a)-(b) . 

Figure 48 (a) -(b) show the Fourier Spectra of Figure 
25 46 (a) - (b) . 

Figure 49 (a) -(b) show the detailed Hilbert Spectra of 
Figure 46(a)- (b) . 

Figure 5 0 shows the wide banded Fourier Spectrogram of 
Figure 45(a), constructed from 128 sampled points. 
30 Figure 51 (a) -(b) show the detailed Fourier Spectrogram of 

Spectra of Figure 46 (a) -(b). 

Figure 52 (a) -(b) show the detailed Morlet Wavelet Spectra 
of Figure 46 (a) -(b). 
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Figure 53 (a) -(b) show the detailed Hilbert Spectra of 
Figure 49 (a) - (b) . 

Figure 54 (a) -(b) show the detailed Hilbert Spectra and 
the acoustical signal data of Figure 46 (a) -(b). 
5 Figure 5 5 shows the acoustical signal data of the sounds, 

"Halloo + Ding," "Halloo, " and "Ding," respectively. 

Figure 56(a) shows the EMD decomposed IMF components of 
sound "Halloo + Ding" , wherein each of the components is 
plotted in the uniform vertical scale. 
10 Figure 56 (b) shows the EMD decomposed IMF components of 

h- sound "Halloo + Ding", wherein each of component is plotted in 
p the vertical scale normalized within its own frame, 
j"*! Figure 57 (a) shows the filtered IMF components of the 

it] sound "Halloo + Ding", wherein the signal associated with the 

= " • 

lis sound "Ding" has been eliminated. 

s Figure 57 (b) shows the IMF components of the sound 

Q 

ry "Ding," separated from Figure 56(a). 

f = ! Figure 58(a) shows the Hilbert Spectrum of the sound 

□ "Halloo + Ding. " 

l % Figure 58(b) shows the HHT filtered Hilbert Spectrum of 

the sound "Halloo + Ding." 

Figure 59(a) shows the acoustical signal data, filtered 
with HHT and Fourier, respectively. 

Figure 59 (a) shows Fourier Spectra for the sound "Halloo 
25 + Ding" and various filtered data. 

Figure 60 (a) -(c) show the acoustical data from a grinder 
operating on a hard surface continuously for 95, 12 0, and 200 
hours , respectively . 

Figure 61 (a) -(c) show the EMD decomposed IMF components 
30 of Figure 60 (a) -(c). 

Figure 62 (a) -(c) show the Hilbert Spectra of Figure 
60 (a) - (c) . 
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DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 

Before describing the computer implemented Empirical Mode 
Decomposition method in detail and its application to 
biological data and acoustical data, the definition and 
physical meaning of intrinsic mode functions in general will 
be discussed. 

Intrinsic Mode Function 

An Intrinsic Mode Function (IMF) is a function that 
satisfies the following two conditions: 

(a) in the whole data set, the number of extrema and 
the number of zero-crossings must either be equal or 
differ at most by one, and 

(b) at any point, the mean value of upper envelope 
defined by the maxima and the lower envelope defined 
by the minima is zero. 

The first condition shares some similarity to the 
traditional narrow band requirements for a stationary Gaussian 
process. The second condition is a totally new idea. 
Conceptually, the second condition modifies the classical 
global requirement to a local one. Furthermore, the second 
condition has the desirable result that the instantaneous 
frequency will not have unwanted fluctuations induced by 
asymmetric wave forms. Mathematically, the second condition 
should ideally be 'the local mean of the data being zero. ' 
For nonstationary data, the x local mean' requires a 'local 
time scale' to compute the mean, which is impossible to 
define. Fortunately, the local time scale need not be defined 
to fulfil the second condition, as will be discussed below. 

To apply these concepts to physical data, the invention 
utilizes the local mean of the signal envelopes to force the 
local symmetry. 

The signal envelopes are defined by the local maxima and 
the local minima. This is an approximation which avoids the 
definition of a local averaging time scale. 



With the physical approach and the approximation adopted 
here, the inventive method does not always guarantee a perfect 
instantaneous frequency under all conditions. Nevertheless, 
it can be shown that, even under the worst conditions, the 
instantaneous frequency so defined is still consistent with 
the physics of the system being studied and represents the 
system being studied much more accurately than previous 
techniques based on Fourier analysis . 

The term "Intrinsic Mode Function" is adopted because it 
represents the oscillation mode embedded in the data. With 
this definition, the IMF in each cycle, defined by the zero- 
crossings, involves only one mode of oscillation. In other 
words, each IMF represents only one group of oscillation modes 
or time scales and no riding waves are allowed. 

Before presenting the inventive EMD method for 
decomposing the data into IMFs, a qualitative assessment of 
the intrinsic oscillatory modes may be roughly determined by 
simply examining the data by eye. From this examination, one 
can immediately identify the different scales directly in two 
ways: the time lapse between the successive alternations of 
local maxima and minima and the time lapse between the 
successive zero-crossings reveals the different scales. The 
interlaced local extrema and zero-crossings give us 
complicated data: one undulation is riding on top of another, 
and they, in turn, are riding on still other undulations, and 
so on. Each of these undulations defines a characteristic 
scale or oscillation mode that is intrinsic to the data: 
hence, the term "Intrinsic Mode Function" is adopted. 

To reduce the data into the needed IMFs, the invention 
utilizes a computer implemented Empirical Mode Decomposition 
Method which is described below. 

Empirical Mode Decomposition (EMD) : The Sifting Process 

First, the Empirical Mode Decomposition method which 
deals with both nonstationary and nonlinear data will be 



discussed. Then, the physical meaning of this decomposition 
will be presented. 

The essence of the EMD method is to identify empirically 
the intrinsic oscillatory modes by their characteristic time 
5 scales in the data, and then decompose the data accordingly. 
The decomposition is based on the following assumptions: 

a. the signal has at least two extrema: one maximum 
and one minimum, and 

b. the characteristic time scale is defined by the 
10 time lapse between the extrema. 

p In other words, the invention adopts the time lapse 

between successive extrema as the definition of the time scale 
U| for the intrinsic oscillatory mode because it gives a much 
Jit finer resolution of the oscillatory modes and because it can 

Ml 

Sis be applied to data with non-zero mean (either all positive or 
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all negative values, without zero-crossings) . A systematic 
way to extract the intrinsic mode functions is the computer 
id implemented Empirical Mode Decomposition method or Sifting 
*zf. Process which is described as follows . 
20 Figure 1(a) illustrates the overall inventive method 

including the Sifting Process in step 120. First, the 
physical activity, process or phenomenon is sensed by an 
appropriate sensor in step 100. Appropriate sensors for 
detecting the physical activity and generating a physical 
25 signal representative thereof are discussed in the practical 
examples below. As an equivalent alternative, the physical 
signal can be inputted in step 100. 

After sensing in step 100, the analog signal is converted 
to the digital domain suitable for computer processing in the 
30 A/D conversion step 105. Depending upon whether the input 
signal is analog or digital step 105 may be bypassed. 

Next, an optional smoothing step 110 may be applied to 
the physical signal. The optional smoothing step 110 may be 



27 



applied to smooth the signal with, for example, a weighted 
running average to remove excessive noise. 

Thereafter, the Sifting Process is applied in step 12 0 to 
Sift the signal with the Empirical Mode Decomposition method 
5 and thereby extract the intrinsic mode function (s). The 
intrinsic mode functions can then be displayed as shown in 
step 130 and checked for orthogonality in step 135. 

Before continuing with the main flow in Figure 1(a), the 
details of the Sifting Process will be explained with 
10 reference to the high level flowchart in Figures 2(a), 2(b) 
12 and the series of graphs showing illustrative results of the 
□ Sifting Process in Figures 3 (a) -(f). 

h l As shown in Figure 1(b), the digitized physical signal 

tf] from step 105 is first windowed by framing the end points in 

in 

gi5 step 107. Then, the Sifting Process begins at step 200 by 

: = . identifying local maximum values of the digitized, framed 

rlJ physical signal from step 107. Figure 3(a) shows a typical 

\~] physical signal 10 which, in this example, represents wind 

Li 5 

O speed spanning a time interval of one second. 

Hi 

20 Before construction of the envelope in steps 210 and 23 0, 

optional intermittency tests (201,221) may be introduced to 
alleviate the alias associated with intermittence in the data 
that can cause mode mixing. 

Optional intermittency test 201 checks the distance 

25 between successive maxima to see if this distance between is 
within a pre-assigned value n times the shortest distance 
between waves. If no, then an intermittency exists and the 
method proceeds to step 202. If yes, then there is no 
intermittency and the upper envelope is constructed in step 

30 210 as further described below. 

Similarly optional intermittency test 221 checks the 
distance between successive minima to see if this distance is 
within a pre-assigned value n times the shortest distance 
between waves. If no, then an intermittency exists and the 
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method proceeds to step 222. If yes, then there is no 
intermittency and the upper envelope is constructed in step 
23 0 as further described below. 

An example of such intermittency is given in Figure 3(g), 
5 in which small scale waves appear only intermittently. By 
strict application of the Sifting Process, the resulting IMFs 
are given in Figures 3(b)-(j), in which two drastically 
different time scales are present in the first IMF component 
as shown in Figure 3 (h) . This mixing of modes breaks up the 
10 main wave train by the small intermittent oscillations. 
Ill If intermittency tests (201,222) are employed which 

□ utilize a preassigned value of n-times the shortest distance 
between waves, the resulting IMFs are shown in Figures 3 (k) - 
?:f (m) , in which the modes are clearly and physically separated. 
-op The effective step to eliminate the mode mixing is achieved by 
~ treating the local extrema which failed the intermittency test 
fll as local maxima and minima (steps 202 and 212), respectively. 

fas a 

Li Therefore, the upper and lower envelope will be identical as 
%J the original data reference line. 

fS = 
I 

20 These intermittency tests (201,221) and the further steps 

(202,222) are optional. By selecting an artificially large n 
value utilized in steps 201 and 221 to test for intermittency, 
the test will be effectively bypassed. Otherwise, the test 
can be bypassed at the initial selection of the program. 

25 Then, the method constructs an upper envelope 2 0 of the 

physical signal 10 in step 210. The upper envelope 20 is 
shown in Figure 3(b) using a dot-dashed line. This upper 
envelope 2 0 is preferably constructed with a cubic spline that 
is fitted to the local maxima. 

30 Next, the local minimum values of the physical signal 10 

are identified in step 220. To complete the envelope, a lower 
envelope 3 0 is constructed from the local minimum values in 
step 230. The lower envelope 30 is also shown in Figure 3(b) 
using a dot-dash line. Like the upper envelope 20, the lower 
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envelope 3 0 is preferably constructed with a cubic spline that 
is fitted to the local minima. 

The upper and lower envelopes 20, 3 0 should encompass all 
the data within the physical signal 10. From the upper and 
5 lower envelopes 20, 30, an envelope mean 40 is the determined 
in step 240. The envelope mean 40 is the mean value of the 
upper and lower envelopes 20, 30. As can be seen in Figure 
3 (b) , the envelope mean 40 bisects the physical signal 10 
quite well. 

io Then, the method generates the first component signal h ± 

I" in step 250 by subtracting the envelope mean 40 from the 



O physical signal 10. This computer implemented step may also 
j.l be expressed as : 



ills Where the envelope mean 40 is m ± and the physical signal is 



H Figure 3(c) shows the first component signal h x . Ideally, 

! J 

the first component signal ^ should be an IMF, for the 
construction of h ± described above seems to have made h ± 

20 satisfy all the requirements of an IMF. In reality, however, 
a gentle hump that resides on a slope region of the data can 
become an extremum when the reference coordinate is changed 
from the original rectangular coordinate to a curvilinear 
coordinate. For example, the imperfection of the envelopes 

25 20, 30 can be seen by observing the overshoots and undershoots 
at the 4.6 and 4.7 second points in Figure 3(b). 

An example of this amplification can be found in the 
gentle hump between the 4.5 and 4.6 second range in the data 
in Figure 3 (a) . After the first round of sifting, the gentle 

3 0 hump becomes a local maximum at the same time location in the 
first component signal Ji a shown in Figure 3(c) . New extrema 
generated in this way actually recover the proper modes lost 
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X(t) - m 1 



= h ± . . . . 1 



5 



X(t) . 
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in the initial examination. Thus, the Sifting Process 
extracts important information from the signal which may be 
overlooked by traditional techniques. In fact, the Sifting 
Process can recover low amplitude riding waves, which may 
appear as gentle humps in the original signal, with repeated 
sittings . 

Still another complication is that the envelope mean 40 
may be different from the true local mean for nonlinear data. 
Consequently, some asymmetric wave forms can still exist no 
matter how many times the data are sifted. This must be 
accepted because, after all, the inventive method is an 
approximation as discussed before. 

Other than these theoretical difficulties, on the 
practical side, serious problems of the spline fitting can 
occur near the ends, where the cubic spline being fit can have 
large swings. Left by themselves, the end swings can 
eventually propagate inward and corrupt the whole data span 
especially in the low frequency components. A numerical 
method has been devised to eliminate the end effects details 
of which will be given later. Even with these problems, the 
Sifting Process can still extract the essential scales from 
the data . 

The Sifting Process serves two purposes: to eliminate 
riding waves and to make the wave profiles more symmetric. 
Toward these ends, the Sifting Process has to be repeated. 
Because only the first component signal h ± has been generated 
so far, the decision step 260, which tests successive 
component signals to see if they satisfy the definition of an 
IMF, is bypassed during the first iteration. 

Thus, step 2 65 is performed which treats the component 
signal as the physical signal in the next iteration. The next 
iteration is then performed by executing steps 200-250. In 
step 250, the second component signal h 1± is generated by 
subtracting the envelope mean from the physical signal (in 



this iteration, the first component signal h x is treated as the 
physical signal). In more formal terms: 

- 2. = ^11 ^ 

Figure 3(d) shows the second component signal h lx . 
Although the second sifting shows great improvement in the 
signal with respect to the first sifting, there is still a 
local maximum below the zero line. After a third sifting, the 
result (third component signal h 12 ) is shown in Figure 3 (d) . 
Now all the local maxima are positive, and all the local 
minima are negative, but to ensure this configuration is 
stable, the Sifting Process should be further repeated. In 
general, the Sifting Process is repeated at least 3 more times 
and, in general, K times to produce h lk . If no more new 
extrema are generated, then h lk is an IMF. In formal terms: 

h l(k-l) ~ m lk = h lk '" 3 

The resulting first IMF component is shown in Figure 3(f) 
after 9 sif tings. The first IMF component of the physical 
signal may be designated as such in step 270 and stored in 
step 275 in memory 415: 

c l = h lk 4 

As mentioned above, all these manipulations are carried 
out numerically in a computer 410. There is not explicit 
close form analytic expression for any of the computer 
implemented steps. 

As described above, the process is indeed like sifting of 
the data by the computer 410 because it separates the finest 
(shortest time scale) local mode from the data first based 
only on the characteristic time scale. The Sifting Process, 
however, has two effects: 
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a. to eliminate riding waves, and 

b. to ensure the envelopes generated by maxima and 
minima are symmetric . 

While the first condition is necessary for the 
5 instantaneous frequency to be meaningful (as discussed below) , 
the second condition is also necessary in case the neighboring 
wave amplitudes have too large a disparity. 

Unfortunately, the effects of the second condition, when 
carried to the extreme, could obliterate the physically 
10 meaningful amplitude fluctuations. Therefore, the Sifting 

Process should be applied with care, for carrying the process 
to an extreme could make the resulting IMF a pure frequency 
modulated signal of constant amplitude. 
*B To guarantee that the IMF component retains enough 

in 

gj5 physical sense of both amplitude and frequency modulations, a 
stopping criterion is employed to stop the generation of the 
next IMF component . 

This stopping criterion is part of the computer 

implemented method and is shown as step 260 in Figure 1(c). 

FU 

20 Step 260 is a decision step that decides whether the stopping 
criteria has been satisfied. The preferred stopping criteria 
determines whether three successive component signals satisfy 
the definition of IMF. If three successive component signals 
all satisfy the definition of IMF, then the Sifting Process 

25 has arrived at an IMF and should be stopped by proceeding to 
step 270. If not, step 260 starts another iteration by 
proceeding to step 2 65 as described above. 

Alternatively, another stopping criteria could be used 
that determines whether successive component signals are 

30 substantially equal. If successive component signals are 

substantially equal, then the Sifting Process has arrived at 
an IMF and should be stopped by proceeding to step 270. If 
not, step 2 60 starts another iteration by proceeding to step 
2 65 as described above. 
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Determining whether successive component signals are 
substantially equal in the alternative stopping criteria 
limits the size of the standard deviation, sd, computed from 
the two consecutive sifting results as follows: 

.5 

5 A very rigorous and preferred value for sd is set between 

0.2 and 0.3. Of course, if faster processing is desired, then 

a trade-off such as a less rigorous value for sd may be used. 

□ Overall, the first IMF component c x should contain the 

§1 finest scale or the shortest period component of the physical 

Mi!o signal 10 . 
v?§ 

\l Before extracting the next IMF component, a test should 

be conducted to determine if the Sifting Process should stop. 

j i 

flj The stopping criteria is shown in Step 3 00. Step 3 00 

-.- s determines whether the component signal has more than 2 

Site extrema. If not, all of the IMF's have been extracted and the 

5 — 

Sifting Process is stopped by proceeding to step 310. If so, 
then additional IMF's may be extracted by continuing the 
process in step 320. 

Step 270 recognizes that an IMF component has been 
20 successfully generated by the Sifting Process by setting the 
component signal equal to an intrinsic mode function. More 
specifically, step 270 causes the computer 410 to store the 
component signal hlk as an intrinsic mode function in memory 
415 . 

25 Then, the first IMF is separated from the physical signal 

in step 290 to generate a residual signal. In particular, a 
residual signal is generated by subtracting the intrinsic mode 
function from the physical signal. In formal terms: 

X(t) - c 1 = r ± .6 



T 

sd = & 
fc=0 
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Because the residue, r lf still includes information of 
longer period components, it is treated as the new physical 
data and subjected to the same Sifting Process as described 
above. Step 32 0 performs this function by treating the 
5 residual signal as the physical signal in the next iteration. 
Thereafter, the next iteration is performed beginning with the 
execution of step 200 as described above. 

The Sifting Process is analogous to a mechanical sieve, 
except it is implemented here in specially programmed computer 
10 and applied to any digital data numerically rather than 
L = v mechanically. 

1% The Sifting Process is repeated for all the subsequent 

"H Tj's. This iterative procedure may be expressed as: 

UJ 

KS. 

m 

S| r i - c 2 = r 2 • 

a 000 ' 7 

f|| r n-l c n ~ r n • 

s . 

iii 

===r 
fll 
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Step 300 stops the Sifting Process by proceeding to stop 
step 310 if the residual signal r D has more than 2 extrema . 
Otherwise, the method proceeds to step 320. 

In other words, Step 310 stops the Sifting Process if the 

20 residual signal r n is monotonically increasing or decreasing. 
This stopping criterion is based on the fact that an IMF 
cannot be extracted from a monotonic function. If r n is not 
monotonically increasing/decreasing, then a next iteration is 
performed beginning with step 320. 

25 Even for data with zero mean, the final residue still can 

be different from zero. For data with a trend, the final 
residue should be that trend. 
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In summary, the Sifting Process decomposes the physical 
signal X(t) into a series of intrinsic mode functions and a 
residue which may be expressed as: 

n 

X(t) = & c 7 < + r n . 8 
i = l XJ 

In other words, the invention extracts a series of IMFs 
5 by Sifting the physical signal with a computer implemented 
Ll Empirical Mode Decomposition method. This IMF series cannot 
□ be generated or derived by any analytic method. It can only 
sa be extracted by the invention through a specially programmed 
computer through the Sifting Process. 

Lio 

EMD Filtering and Curve Fitting 

Figure 1 (d) illustrates the inventive method of 

n i 

11 performing filtering and curve fitting. The method begins, as 

Yd in Figure 1(a), with sensing the physical phenomenon or 

figs inputting a signal representative thereof in step 100. The 

optional smoothing 105 may then be applied to eliminate excess 

noise . 

Then, the signal is Sifted with the Empirical Mode 
Decomposition Process to extract the intrinsic mode functions 

20 in step 120. The IMFs may then be directly displayed 130, 
stored 132 or transmitted 134. 

As mentioned above, fitting a curve to a signal is not 
always possible. The reasons include too many data points for 
efficient curve fitting and nonlinear data that does not admit 

25 an acceptable curve fit. The IMFs, however, extract important 
information from the original signal while reducing the number 
of data points and simplifying the resulting signal. By 
properly selecting only the most relevant IMFs, one can 
construct a filtered version of the original signal. 
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Step 350 performs this process of constructing a filtered 
signal from a subset of the IMFs. Selecting which IMFs should 
in the subset and which should not may be an empirical or 
intuitive process. Typically, the IMFs having the lowest 
5 frequency components are chosen as part of the selected subset 
with the higher frequency components excluded therefrom. In 
this way, a simplified (filtered) version of the original 
signal can be constructed which eliminates unnecessary and, 
hopefully, not physically components (IMFs). 
10 More specifically, step 350 performs a summation of the 

M selected IMFs. The selection of IMFs can be automated, e.g. 
p\ the n lowest frequency IMFs, or manual. 

s ^ Generating the filtered signal typically permits curve 

LiJ 

fitting to be successful and efficient. Once the filtered 
\Hs signal is generated, conventional curve fitting techniques 
a (step 360) such as a least squares estimation process can be 
fll utilized. 

H The fitted curve may then be outputted 3 62, displayed 

rs 3 64, stored 3 66 or transmitted 3 68. 

T= = S- 

as : 

r 5'0 The filtered signal may also be subjected to further 

processing in step 370. Namely, the oscillatory energy about 
the mean is calculated in step 370. This calculation is 
essentially the same as the instantaneous energy density 
calculation explained below, with the calculation being 

25 performed for a subset (k) of the IMFs. Further details and 
an example are discussed below. The result of step 37 0 is a 
signal that may be subjected to the Sifting Process to extract 
a set of IMFs in step 120. Although Figure 1(d) illustrates 
an infinite loop (steps 120, 350 and 370), at some point it 

30 will not be appropriate or possible extract further IMFs. 
Thus, this loop is not really infinite and is typically 
performed one or two times . 

Computer for Implementing Inventive Method 
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A computer suitable for programming with the inventive 
method is diagrammat ically shown in the block diagram of 
Figure 2. The computer 410 is preferably part of a computer 
system 400. To allow human interaction with the computer 410, 
5 the computer system includes a keyboard 430 and mouse 435. 

The computer programmed with the inventive method is analogous 
to a mechanical sieve: it separates digital data into series 
of IMF's according to their time scales in a manner analogous 
to a mechanical sieve which separates aggregated sand 
10 particles according to their physical size. 

I* Because the invention is applied to analyze physical 

□ signals, the computer system 400 also includes an input device 
440, sensor 442 and/or probe 444 which are used to sample a 

*f) physical phenomenon and generate physical signal 

ill 

^5 representative thereof. More particular examples of such 

= =i inputs (440, 442 and 444) are described in relation to Figures 

Ki 21-25. 

H To output the results of the computer implemented method, 

iU 

O the computer system 400 also includes a display 450 such as a 

til 

*5o cathode ray tube or flat panel display, printer 460 and output 
device 470. Each of these outputs (450, 460, 470) should have 
the capability to generate or otherwise handle color outputs 
because, for example, the Hilbert Spectrum may be in color. 
The generalized output device 470 may also include a 

25 network connection to connect the computer 400 to a local or 
wide area network. In this way, the physical signal may be 
inputted from the network. Furthermore, all outputs can be 
sent to another location via such a network connection. 

Furthermore, the computer system 40 0 also includes a mass 

30 storage device 420. The mass storage device 420 may be a hard 
disk, floppy disc, optical disc, etc. The mass storage device 
42 0 may be used to store a computer program which performs the 
invention when loaded into the computer 410 . As an 
alternative, the input device 440 may be a network connection 
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or off-line storage which supplies the computer program to the 
computer 410. 

More particularly, the computer program embodiment of the 
invention may be loaded from the mass storage device 42 0 into 
5 the internal memory 415 of the computer 410. The result is 
that the general purpose computer 410 is transformed into a 
special purpose machine that implements the invention. 

Even more particularly, each step of inventive method 
will transform at least a portion of the general purpose 
10 computer 410 into a special purpose computer module 
H implementing that step. For example, when the sifting step 

a 

q 12 0 is implemented on the computer 410, the result is a 

y m computer implemented sifting apparatus (sifter) that performs 

W 

% Q the sifting functions of sifting step 120. 

ris Other embodiments of the invention include firmware 

* embodiments and hardware embodiments wherein the inventive 

i\l method is programmed into firmware (such as EPROM, PROM or 

3 %T 

PLA) or wholly constructed with hardware components. 

Ill 

p Constructing such firmware and hardware embodiments of the 
s 3?o invention would be a routine matter to one of ordinary skill 
using known techniques. 

Article of Manufacture 

Still further, the invention disclosed herein may take 
25 the form of an article of manufacture. More specifically, the 
article of manufacture is a computer-usable medium, including 
a computer-readable program code embodied therein wherein the 
computer-readable code causes computer 410 to execute the 
inventive method. 
30 A computer diskette such as disc 480 in Figure 2 is an 

example of such a computer-usable medium. When the disc 480 
is loaded into the mass storage device 480, the computer- 
readable program code stored therein is transferred into the 
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computer 410. In this way, the computer 410 may be instructed 
to perform the inventive methods disclosed herein. 

Illustration of Sifting Process 

5 To further illustrate the Sifting Process, consider 

Figure 4(a) which shows a physical signal representing wind 
data collected in a laboratory wind-wave tunnel with a high 
frequency response Pitot tube located 10 cm above the mean 
water level. The wind speed was recorded under the condition 

10 of an initial onset of water waves from a calm surface. 

ilj Clearly, the physical signal is quite complicated with many 

Q local extrema but no zero-crossings such that the time series 

SI 

m represents all positive numbers. 

)li Although the mean can be treated as a zero reference, 

LI 5 

-aj5 defining it is difficult, for the whole process is transient. 
This example illustrates the advantage of adopting the 

is! 

ill successive extrema for defining the time scale and the 

s : 
?« 

jjj difficulties of dealing with nonstationary data. In fact, a 
id physically meaningful mean for such data is impossible to 

p § 
3 - - 

20 define using standard methods. The EMD eliminates this 
difficulty. 

Figure 4(b) shows the wind speed signal of Figure 4(a) on 
a different scale for comparison purposes. Figures 4(c)-(k) 
show all the IMFs obtained from repeatedly sifting the wind 

25 speed signal in Figure 4 (b) . The efficiency of the invention 
is also apparent: the Sifting Process produces a total of 9 
intrinsic mode function components while the Fourier transform 
needs components which number as many as half of the total 
number of points in a given window to represent the wind data 

3 0 with the same accuracy. 

The separation of the wind speed data into locally non- 
overlapping time scale components is clear from Figures 4(c)- 
(k) . In some components, such as c 2 and c 3 , the signals are 
intermittent, then the neighboring components might include 
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oscillations of the same scale, but signals of the same time 
scale would never occur at the same locations in two different 
IMF components . 

The components of the EMD are usually physical, for the 
5 characteristic scales are physically meaningful. 

To confirm the validity and completeness of the Sifting 
Process, the wind speed signal can be reconstructed from the 
IMF components. Figures 5(a)- (i) show this reconstruction 
process starting from the longest period IMF to the shortest 
10 period IMF in sequence. For example, Figure 5(a) shows the 
M wind speed signal and the longest period component, c 9 , which 
q is actually the residue trend, not an IMF. 

By itself, the fitting of the trend is quite impressive, 
%li and it is very physical: the gradual decrease of the mean wind 

; ==? 

*3|5 speed indicates the lack of drag from the calm surface 
= s initially and the increasing of drag after the generation of 
m wind waves. As the mean wind speed deceases, the amplitude of 
the fluctuation increases, another indication of wind-wave 

U.J 

Q interactions. 

By adding the next longest period component, c 8 , the trend 
of the sum, c 9 + c 8 , takes a remarkable turn, and the fitting 
to the wind speed signal looks greatly improved as shown in 
Figure 5 (b) . Successively adding more components with 
increasing frequency results in the series of Figures 5(c)- 

25 (i). The gradual change from the monotonic trend to the final 
reconstruction is illustrative by itself. By the time the sum 
of IMF components reaches c 3 in Figure 5 (g) , essentially all 
the energy contained in the wind speed signal is recovered. 
The components with the highest frequencies add little more 

30 energy, but they make the data look more complicated. In 
fact, the highest frequency component is probably not 
physical, for the digitizing rate of the Pitot tube is too 
slow to capture the high frequency variations. As a result, 
the data are jagged artificially by the digitizing steps at 
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this frequency. The difference between the original data and 

the re-constituted set from the IMF 1 s is given in Figure 5(j). 

The magnitude of the difference is lCT 15 , which is the limit of 
the computer 410. 



Having obtained the Intrinsic Mode Function components 
(through either the local extrema or curvature extrema Sifting 
Processes) , the next step in the computer implemented method 
is to apply the Hilbert Transform to each component and 
generate the Hilbert Spectrum as shown in step 140 in Figure 
1(a) . A brief tutorial on the Hilbert transform with 

emphasis on its physical interpretation can be found in Bendat 
and Piersol, 1986, "Random Data: Analysis and Measurement 
Procedures", 2nd Ed., John Wiley & Sons, New York, NY. 
Essentially, the Hilbert transform is the convolution of X(t) 
with 1/t. This convolution emphasizes the local properties of 
X(t) . In more formal terms, given a time series X(t) , the 
Hilbert Transform Y(t) can be expressed as 



where P indicates the Cauchy principal value. 

With this definition, X(t) and Y(t) form a complex 
conjugate pair. This complex conjugate pair Z(t) may be 
expressed as: 



The Hilbert Spectrum 




— oo 



Z(t) 



X(t) + iY('t) = a(t) e i( 3 (t) ,10 



in which 



1 

a(t) = [ X 2 (t) + Y 2 (t) ] 2,11 



q(t) = arctanyj t ) . 



X(t) 



12 
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After performing the Hilbert transform to each IMF 
component, we can express the time series data X(t) in the 
following form: 

n 

X(t) = & a..(t) U Wj'(t) dt 13 

j = l ^ 

In Equation (13), the residue, r nf is purposefully 
omitted, for it is either a monotonic function, or a constant. 
Although the Hilbert transform can treat the monotonic trend 
as part of a longer oscillation, the energy involved in the 
residual trend could be overpowering. In consideration of the 
uncertainty of the longer trend, and in the interest of 
information contained in the other low energy and higher 
frequency components, the final non-IMF component should be 
left out. It, however, could be included, if physical 
considerations justify its inclusion. 

Note that Equation (13) gives both amplitude and 
frequency of each component as functions of time. It should 
be pointed out that no analytical method can generate the 
expression in Equation (13). Instead, all the components may 
be extracted only by a specially programmed computer applying 
the inventive Sifting Process and the Hilbert transform. The 
variable amplitude and frequency have not only greatly 
improved the efficiency of the expansion, but also enabled the 
expansion to accommodate nonstationary data. With IMF 
expansion, the amplitude and the frequency modulations are 
also clearly separated. 

Equation (13) also enables the computer implemented 
method to represent the amplitude and frequency as functions 
of time in a three-dimensional plot, in which the amplitude 
can be contoured on the frequency- time plane. This frequency- 
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time distribution of the amplitude is designated as the 
Hilbert Amplitude Spectrum, H(_, t) , or simply Hilbert 
Spectrum. Thus we have: 

n 

JT(w, t) = A a 7 '(t)e i ^w :f (t)dt 14 

j = i J 

In which Hf_, t>) is the Hilbert spectrum of the frequency (_) 
and time (t) and a^t) is the j-th component of the IMF. 

In the presentation, the amplitude (with or without 
smoothing) can be expressed in color maps, black-grey maps, or 
contour maps. Color maps, however, greatly increase the 
operator's ability to fully analyze the spectrum. In some 
cases, a color map will permit the operator to discern 
relationships and trends that would not be apparent in black- 
grey maps thereby making a color display a necessary component 
in some cases. 

If amplitude squared is more desirable to represent 
energy density, then the squared values of amplitude can be 
substituted to produce a Hilbert Energy Spectrum just as well. 

Various forms of Hilbert spectra presentations can be 
generated by the computer in the display step 190: both color 
coded maps and contour maps may be employed to present the 
Hilbert spectra with or without smoothing. The Hilbert 
Spectrum in the color map format for the wind data is shown in 
Figure 6(a). The Hilbert spectrum in Figure 6(a) gives a 
very different appearance when compared with the corresponding 
Wavelet spectrum shown in Figure 6(b). While the Hilbert 
Spectrum in Figure 6(a) appears only in the skeleton form with 
emphasis on the frequency variations of each IMF, the Wavelet 
analysis result gives a smoothed energy contour map with a 
rich distribution of higher harmonics. 

If a more continuous form of the Hilbert Spectrum is 
preferred, a smoothing method can be optionally applied in 
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step 155. The first type of a smoothing method which may be 
used in the invention is a weighted spatial filter which 
averages over a range of cells. For example, a 15 by 15 
weighted Gaussian filter may be employed in step 155 as is 
5 known in the art to smooth this data. Figure 6(c) shows the 
result of applying the 15 by 15 weighted Gaussian filter. 

Although smoothing step 155 degrades both frequency and 
time resolutions, the energy density and its trends of 
evolution as functions of frequency and time are easier to 
10 identify. In general, if more quantitative results are 
M desired, the original skeleton presentation is better. If 

a 

p more qualitative results are desired, the smoothed 

presentation is better. As a guide, the first look of the 

ty 

»Q data is better in the smoothed format. 

i n 

IS|5 The alternative of the spatial smoothing is to select a 

5 lower frequency resolution and leave the time axis 

Q 

ry undisturbed. The advantages of this approach are the 
l* preservation of events' locations and a more continuous 
ill frequency variation. Furthermore, a lower frequency 
'20 resolution saves computation time for the computer implemented 
method. 

To optimize such computation time, the optimal frequency 
resolution in the Hilbert spectrum can be computed as follows. 
Let the total data length be T, and the digitizing rate of the 

25 sensor be _t. Then, the lowest frequency that can be 

extracted from the data is 2/T Hz, which is also the limit of 
frequency resolution for the data. The highest frequency that 
can be extracted from the data is l/(n _t) Hz, in which n 
represents the minimum number of _t needed to define the 

30 frequency accurately. 

Because the Hilbert transform defines instantaneous 
frequency by differentiation, more data points are needed to 
define an oscillation. The absolute minimum number of data 
points is five for a whole sine wave. Although a whole sine 
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wave is not needed to define its frequency, many points within 
any part of the wave are needed to find a stable derivative. 
Therefore, the maximum number of the frequency cells, W, of 
the Hilbert spectrum should be 

1 

nDt T 
N = ~T~ = rTDt " 
T 

In order to make the derivative stable, the data is 
averaged over three adjacent cell values for the final 
presentation . 

To illustrate, consider the wind data of Figure 4(a) 
which was digitized at a rate of 0.01 seconds and has a total 
length of 30 seconds. Therefore, the highest frequency that 
can be extracted is 25 Hz. The total cell size could be 600, 
but they have been averaged to 200 in Figure 6(a) . 

Marginal Spectrum 

The marginal spectrum offers a measure of total amplitude 
(or energy) contribution from each frequency value. In other 
words, the marginal spectrum represents the cumulated 
amplitude over the entire data span. 

As shown in Figure 1(a), the marginal spectrum is 
calculated by the computer implemented method in step 145 
after applying the Hilbert Transform in step 140. The 
marginal spectrum is the Hilbert Spectrum integrated through 
all time. In this simplification, the time coordinate is lost 
as in the Fourier spectral analysis, which leaves a summary of 
the frequency content of the event. Therefore, this 
presentation should only be used when the phenomena being 
analyzed is stationary. Formally, the marginal spectrum h(_) 
is defined as: 
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T 

h (w) = U h (w, t ) dt .16 
0 




Because there is no analytic expression for H(_,t), the 
integration can only be performed in a computer as a sum. 

An example of a marginal spectrum is shown in Figure 7. 
5 More particularly, Figure 7 shows the corresponding marginal 
spectrum of the Hilbert Spectrum given in Figure 6(a) . 

The frequency in either H(_ r t) or h(_) has a totally 
^ different meaning from results generated by applying Fourier 
W spectral analysis. In the Fourier representation, the 
yyp existence of energy at a frequency, _, means a component of a 

:zi sine or a cosine wave persisted through the time span of the 

y § 

%j data . 

f=! In contrast, the existence of energy at the frequency, 

! y means only that, in the whole time span of the data, there is 

lip a higher likelihood for such a wave to have appeared locally. 

Z:i In fact, the Hilbert Spectrum is a weighted, non-normalized, 

in 

joint amplitude-frequency- time distribution. The weight 
assigned to each time- frequency cell is the local amplitude. 
Consequently, the frequency in the marginal spectrum indicates 

20 only the likelihood of an oscillation with such a frequency 
exists . The exact occurrence time of that oscillation is 
given in the full Hilbert spectrum. 

To illustrate this difference, the corresponding Fourier 
Spectrum of the wind speed signal is also given in Figure 7 

25 using a dotted line. As can be observed, there is little 
similarity between the Fourier spectrum and the marginal 
spectrum. While the Fourier spectrum is dominated by the DC 
term because of the non-zero mean wind speed, the marginal 
spectrum gives a nearly continuous distribution of energy. 

30 The Fourier spectrum is meaningless physically, because the 
data is not stationary. In contrast, the marginal spectrum 
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provides a physically meaningful interpretation of the wind 
speed signal. 

Instantaneous Frequency 

5 There are two types of frequency modulation: the inter- 

wave and the intra-wave modulations. The first type is 
familiar because the frequency of the oscillation gradually 
changes as the waves disperse. Technically, in dispersive 
waves, the frequency is also changing within one wave, but 
10 that is generally not emphasized either for convenience, or 

f h for lack of a more precise frequency definition. The second 

P 

El type is less familiar, but it is also a common phenomenon: if 

the frequency changes from time to time within a wave its 
'%Q profile can no longer be a simple sine or cosine function. 

in 

ojs Therefore, any wave profile deformation from the simple 

1. sinusoidal form implies intra-wave frequency modulation. 

fy In the past, such phenomena were treated as harmonic 

\~] distortions. Nevertheless, such deformations should be viewed 

Hi as intra-wave frequency modulation because the intra-wave 

'20 frequency modulation is a more physically meaningful term. 

In order to understand these frequency modulations, the 
invention applies a unique definition of instantaneous 
frequency. This definition stems from the EMD method and 
requires the signal to be reduced into IMF components. After 
25 extracting the 

IMF components, an instantaneous frequency value can be 
assigned to each IMF component. Consequently, for complicated 
data in which more than one IMF may be extracted by EMD, there 
can be more than one instantaneous frequency at a time 
30 locally. 

With the Hilbert Transform, a unique definition of 
instantaneous frequency may be applied by the computer 
implemented method as illustrated by step 160. Step 160 
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calculates the instantaneous frequency which is formally 
defined as follows : 

dq( t) 
w < fc > = "dF - - 17 

By calculating instantaneous frequency, step 160 of the 
invention permits the frequency value to be defined for every 
point with the value changing from point to point. 

The validity and the implications of the instantaneous 
frequency for nonlinear signals may be analyzed by examining 
an idealized Stokes wave in deep water. The wave profile of 
such a Stokeian wave is modeled by 

X( t) =cos (wt+esinwt) 18 

Therefore, it is a intra-wave frequency modulated signal. 
Approximately, equation (18) can be shown to be: 

X( t) = (1+e/ ) cosw t+ecos2wt+ 19 

The wave profile is also shown in Figure 9 (a) . Because the 
intra-wave frequency can only be approximated by harmonics in 
Fourier analysis, we can still have the same profile, but not 
the same frequency content . The wave form shows sharpened 
crests and rounded off troughs, which make the crests and 
troughs asymmetric with respect to the mean surface. 

Processed with computer implemented EMD, this data yields 
only one IMF component as shown in Figure 9(b), with a 
constant offset component (not shown) . Although this wave 
has only one characteristic scale or IMF, the Wavelet analysis 
result shown in Figure 9(c). Figure 9(c) has many harmonics 
with two visible bands of energy corresponding to the highest 
order of approximations of the possible harmonics. 



49 



In contrast, the IMF data can be processed by the 
inventive method to give the Hilbert Spectrum shown in Figure 
9(d) . The Hilbert Spectrum has only one frequency band 
centered around 0.03 Hz, the fundamental frequency of the wave 
train, but there is an intra-wave frequency modulation with a 
magnitude range of 0.02 to 0.04 Hz as the model truly 
represents. This intra-wave frequency modulation has been 
totally ignored in the past, for the traditional definition of 
frequency is based on the reciprocal of periodicity and 
Fourier Analysis. 

Instantaneous Energy Density 

Furthermore, the computer implemented method may also 
calculate the instantaneous energy density in step 150. The 
instantaneous energy density, IE, is defined as: 

IE( t) = U w H 2 (w, t ) dw 20 

Still further, this instantaneous energy density depends 
on time. Therefore, the IE can be used to check energy 
fluctuations . 

Stationarity 

To quantitatively measure the stationarity of a physical 
signal, the invention utilizes step 165 to calculate various 
stationarity measurements. Before introducing the preferred 
stationarity measurements, a brief review of conventional 
stationarity measurements is presented. 

The classic definitions of stationarity are dichotomous : 
a process is either stationary or nonstationary . This crude 
definition is only qualitative in nature. Such definitions 
are both overly stringent and useless at the same time: few 
data sets can satisfy the rigor of these definitions; 
consequently, no one even bothers using them to check 
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stationarity of the signal. As a result, data as 
nonstationary as earthquake and seismological signals are 
routinely treated as stationary (see, for example, Hu, et al . , 
1996., Earthquake Engineering. Chapman & Hall, London). 

Sometimes, for some obviously nonstationary data, two 
less stringent definitions are invoked: piece-wise stationary 
and asymptotically stationary. These definitions are still 
dichotomous . 

To quantify the statistical processes further, an index 
is needed to give a quantitative measure of how far the 
process deviates from stationarity. A prerequisite for such a 
definition is a method to present the data in the frequency- 
time space. 

With the energy- frequency- time distribution (Hilbert 
Spectrum) described above, stationarity of the physical signal 
may be quantitatively determined. Therefore, the invention 
introduces an index of stationarity as follows and calculates 
a Degree of Stationarity in step 165. 

The first step in defining the Degree of Stationarity, 
DS(_), is to find the mean marginal spectrum, n(_ ), as 

1 

n(w) = J, h(w) 21 



Then, the Degree of Stationarity may be defined as: 



D5 (w) = r U ^ i - n(w) J2 dtf ........ 22 



Again, the value of DS (_) can be determined by the computer. 
Therefore, the specialized computer 410 according to the 
invention can be treated as a stationarity meter. 

For a stationary process, the Hilbert spectrum cannot be a 
function of time. Then, the Hilbert Spectrum will consist of 



only horizontal contour lines and DS(_) will then be 
identically zero. Only under this condition will the marginal 
spectrum be identical to the Fourier spectrum, then the 
Fourier spectrum will also make physical sense. On the other 
5 hand, if the Hilbert spectrum depends on time, the index will 
not be zero, then the Fourier spectrum will cease to make 
physical sense. 

In general, the higher the index value, the more 
nonstationary is the process. The DS for the wind data is 
10 shown in Figure 8 (a) . As the index shows, the data are highly 
nonstationary especially for the high frequency components. 

Eq. (22) defines the stationarity as a function of 
frequency. This is necessary because certain frequency 
components can be nonstationary while other components remain 



in 

%|5 stationary. An example is sporadic riding wind-generated 
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waves on otherwise uniform swell: the low frequency swell 
component is stationary, while the high frequency wind waves 
i are intermittent, and hence nonstationary. 

J The degree of stationarity can also be a function of time 

20 implicitly, for the definition depends on the time length of 
integration in Eq. (22) . Therefore, a process can be piece- 
wise stationary. On the other hand, for a singular outburst 
in an otherwise stationary signal, the process can be regarded 
as almost stationary if a long time integral is performed, but 

25 nonstationary the integral only encompasses the immediate 
neighborhood of the outburst. 

Stationarity can be a complicated property of the 
process : for any signal shorter than a typical long wave 
period, the process may look transient. Yet as the signal 

30 length gets longer, the process can have many longer wave 

periods and becomes stationary. On the other hand, the signal 
can be locally stationary while in a long time sense 
nonstationary. An index is therefore not only useful but also 
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necessary to quantify the process and give a measure of the 
stationarity . 

The invention also calculates a Degree of Statistic 
Stationarity in step 165. The degree of stationarity defined 
5 in Eq. (22) can be modified slightly to include statistically 
stationary signals, for which the Degree of Statistic 
Stationarity, DSS(_,_T), is defined as 

1 TT T ( jf(W, t) \ 

DS(w) = y U ^ i - — — p dt 23 

Q 

r\ where the over-line indicates averaging over a definite but 
S J0 shorter time span, _T, than the overall time duration of the 

i:j 

%y data, T. For periodic motions, the _T can be the period. 

til 

s\ Such a time scale, however, is hard to define precisely for 

= high dimensional, nonstationary dynamic systems. 

O 

f|j Even with this difficulty, the definition for DSS could 

Ft5 be more useful in characterizing random variables from natural 

yj 

□ phenomena. Furthermore, DSS will depend on both frequency and 
the averaging time span. For the wind data taken as an 
example, the DSS is given in Figure 8(a) with _T = 10, 50, 
100, and 3 00 time steps respectively. The results show that 

20 while the high frequency components are nonstationary, they 

can still be statistically stationary. Two frequency bands at 
7 and 17 Hz are highly nonstationary as the DSS averaged at 
100 time steps shown in Figure 8(b). These components are 
intermittent as can be seen in the IMF components and the 

25 marginal spectrum. A section of the original wind data is 
also plotted in Figure 8(c) to demonstrate that there are 
indeed prominent 7 and 17 Hz time scale oscillations. 

Display of Selected Results 

30 The invention displays various results of the above- 

described computer implemented method in step 190. These 

53 



displays are extremely useful in analyzing the underlying 
physics of the physical phenomenon being studied as described 
above. Furthermore, particular examples of these displays and 
the increased understanding of the underlying physics which 
these displays permit are discussed in the following section. 

For example, the invention generates various Hilbert 
spectra displays in the display step 190. As mentioned above, 
both color coded maps and contour maps may be employed to 
display the Hilbert spectra in display step 190. In addition, 
the color coded maps convey information to the operator in a 
uniquely accessible way permitting a more thorough and deeper 
understanding of the physical phenomenon and may be considered 
as necessary to analyze some physical phenomena. 

The displays generated by the invention in display step 
190 are central to the invention because they allow an 
operator to analyze the underlying physics of the physical 
phenomenon being studied. 

The display step 190 outputs displays to display 450. As 
mentioned above, display 450 includes devices such as a 
cathode ray tube and a flat panel display. As an alternative, 
display step 290 may generate a hard copy output by utilizing 
printer 4 60 or may send the generated display to output device 
470. 

BIOLOGICAL SIGNAL PROCESSING 

The above description provides a foundation of signal 
analysis methodologies which has broad applicability to a wide 
variety of signal types . Once such applicable signal type is 
biological signals which are signals measured or otherwise 
representative of a particular biological phenomenon. 

For example, any given physiological process has a number 
of quantities that can be measured to produce a signal or 
otherwise be represented by a signal. The invention can be 
applied to such signals in an effort to better understand the 
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underlying biological phenomenon. For example, the 
relationships between different biological variables can be 
studied with a precision and depth of understanding not 
possible before. 

The invention can also be used to arrive at an analytical 
function representative of the biological phenomenon. By 
using the inventive filtering methods, the biological signal 
can be simplified to eliminate components not particularly 
relevant to the analysis and to thereby make it possible to 
represent the signal with a function. 

Blood pressure data, such as the data shown in Figures 
9 (a) -(g), are examples of the invention's applicability to 
analyzing and characterizing biological signals. This blood 
pressure data and the analysis that follows focusses on 
changes in blood pressure over the course of day. 

To collect this data, a pressure probe was implanted into 
the pulmonary artery of a rat. Measurements were taken with 
the rat under controlled conditions (e.g. rat breathing normal 
atmosphere at sea level, quiet room, lit for 12 of the 24 hour 
time span, etc.). In the experiment, a Statham P23ID 
transducer was utilized, data collected by the computer 400, 
and also recorded on a chart recorder. A standard A/D 
converter was utilized to converted the transduced analog 
signal to digital data for storage and processing by the 
computer 400. It is to be understood that this particular 
set-up is exemplary in nature and illustrates one of the 
various ways in which the biological data can be obtained for 
processing by the invention. 

Figure 9(a) shows the results of measuring the rat's 
blood pressure for a 24-hour period with the pressure waveform 
being sampled at a rate of 100 points/second. A one hour 
strip of this data is shown in Figure 9(b). Two separate 10- 
second strips of this data are shown in Figures 9(c) -(d) with 
9(c) exhibiting less regularity than 9(d). Figures 9(e)-(f) 
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are the systolic peaks and diastolic troughs for the one hour 
record shown in Figure 9 (b) , obtained by connecting the 
successive peaks and successive valleys, respectively. 

In order to compare the invention against conventional 
5 methodologies, Figure 10(a) is presented which shows the 

conventional Fourier Spectrum for the one-hour data of Figure 
9 (b) . Figure 10 (a) shows this Fourier Spectrum from which one 
can easily identify spectral peaks at 1.5, 6.5 and 13 Hz. The 
6.5 Hz peak represents the heartbeats which necessarily affect 
10 blood pressure. The 13 Hz peak is the harmonic of the 6.5 Hz 
Ill heartbeat peak. The 1.5 Hz peak is probably related to 
O respiration. 

Ill Figures 10(b) -(c) show conventional Fourier spectra for 

the two 10-second sections of data (Figures 9(c) -(d)). 

in 

"-lis Fourier analysis for the one-hour data of Figure 9 (b) is shown 
in Figure 10 (d) in which the locations of the highest peaks of 
ril the Fourier spectra in every 1-minute window are plotted on 

|ssb 

hj the time- frequency plane. A perspective view of the windowed 

w 1-hour and 10-sec Fourier results Fourier results are given in 

rn 

20 Figures 10(e) and 10(f), respectively and show the variation 
of the amplitudes of the signals. 

Figure 10(g) directly compares the inventive results 
against conventional techniques. In Figure 10(g), the Fourier 
(dotted line) spectra and the Marginal Hilbert (solid line) 

25 spectra of the blood pressure data from Figure 9(c) are 

plotted together. As can be seen from this comparison, the 
Fourier spectra devotes a great deal of energy to the higher- 
harmonic components . 

As explained above, the invention adopts the spacing of 

30 extrema values to analyze the data. Specifically, the sifting 
process is applied to the data to extract a set of intrinsic 
mode functions. 

Figures 11 (a) -(h) show the results of applying the 
empirical mode decomposition method of the invention to the 
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blood pressure data of Figure 9(c). As shown therein, the 
invention produces eight IMFs from this blood pressure data. 
The residual IMF R 8 is constant. 

Figures 12 (a) - (h) show the results of applying the 
Empirical Mode Decomposition method of the invention to the 
blood pressure data of Figure 9(d). As shown therein, the 
invention produces eight IMFs (C^Cq) from this blood pressure 
data. The residual IMF C 8 (a.k.a R 8 ) is constant. 

In Figures 12 (a) - (h) , the IMF components have very 
different amplitudes. The IMFs having the most energy are C 2 , 
C 3 and C 4 . The amplitude and periodicity of these three main 
components closely approximate their respective, constant 
levels . 

By recognizing the main or most significant IMFs, one can 
reconstruct a filtered version of the signal that eliminates 
undesired components while continuing to faithfully represent 
the original signal. For example, Figure 13(b) shows a 
combination of IMFs C 2 , C 3 and C 4 . The sum of these components 
faithfully represents the original signal as can be seen by 
comparing the IMF sum in Figure 13 (b) with the original signal 
in Figure 9(c). 

Of course, other filters can be constructed with the 
IMFs.. One such example is Figure 13(a) which shows a 
summation of IMFs C lf C 2 and C 3 that were calculated from the 
blood pressure signal of Figure 9(c). 

Another filtering example is shown in Figure 13(c) which 
is a sum of IMFs C 4/ C 5 , C 6 , C 7 and C 8 shown in solid line. 
These components are the lower- frequency components and their 
sum recovers the slow variation of the pressure signal. The 
dotted line in Figure 13(c) shows the original pressure 
signal . 

Figure 14 (a) -(b) show the Hilbert spectrums of the IMF 
components derived from the Figure 9(d) and 9(c) blood 
pressure signals, respectively. Specifically, The Hilbert 



57 



spectrum of Figure 14(a) shows that the most prominent energy 
bands are centered at 6.5, 3.0 and 1 . 5 Hz . As shown, there 
are intra-wave frequency modulations in this spectrum which 
are indications of nonlinear dynamics. The wide fluctuations 
5 of frequency values in the Figure 14 (b) Hilbert spectrum make 
any visual mean estimation exceedingly difficult. 

By further applying the inventive methodologies, a 
marginal Hilbert Spectrum can be calculated by integrating the 
Hilbert Spectrum. An example is shown in Figure 10(g) which 

10 uses a solid line to plot the marginal Hilbert spectrum of the 
Figure 9(c) blood pressure signal. Prominent spectral peaks 

p occur at 7 and 1 . 3Hz . 

IZJ 

sj Figure 10 (g) also uses a dotted line to plot the 

W corresponding Fourier Spectrum. By comparing these spectra, 
Uf one can see that the mean peaks line up. However, the Hilbert 
spectrum clearly depicts the fluctuation of energy over 

□ frequencies without allowing the frequency of oscillation to 

ft I 

11 be variable in the whole time window. 

i . * 

One basic difference between the conventional spectra of 
|$p Figure 10(d) and the inventively derived spectra of Figure 

14(a) lies in the stationarity hypothesis. In Figure 10(d), 
the Fourier spectrum assumes stationary oscillation. In 
Figure 14(a), the Hilbert spectra, does not make this 
assumption and is valid for nons tationary oscillations. 
25 By using the invention, one can more readily and 

accurately recognize the various features of blood pressure 
signals. Moreover, the invention offers a more comprehensive 
view of the blood pressure fluctuation than the classical 
Fourier analysis. 
30 The invention may also be used to study the effect of one 

variable on another variable. Specifically, by changing one 
variable of a system, measuring the effect on another 
variable, and then applying the invention one can arrive at a 
much deeper understanding of underlying system. The invention 
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also permits the modelling or representation of data with an 
analytic function that was not possible with conventional 
techniques . 

Biological systems are examples of systems which can be 
5 better understood and modelled by applying the invention. 

Biological systems are often studied by changing one variable 
and recording the changes of other variables The resulting 
data is often oscillatory, stochastic and non-stationary. 
Data having these properties are particularly susceptible to 
10 processing by the invention. 

t _: A particular example elaborated upon here is the effect 

O of changing the breathing gas concentration of oxygen on the 
%j blood pressure of an animal. Analyzing such stochastic data 
-~ to obtain crystal clear results describing the effects of 
Life oxygen concentration changes on blood pressure has been 
impossible with conventional techniques such as Fourier 
1:1 analysis. 

| == Concrete data illustrating this example is shown in 

1% Figures 15 (a) -(d). Specifically, the pulmonary blood pressure 
FJjb in the arterial trunk of two rats breathing normal air at sea 
level is shown in the upper trace of Figures 15 (a) -(b) with 
the lower trace showing the oxygen concentration. 

These same rats where then subjected to step changes in 
oxygen concentration: Figure 15(c) shows a step increase and 
25 15 (d) shows a step decrease in oxygen concentration and the 

effect thereof on blood pressure. From this raw data, one can 
see certain trends. However, there is no definitive way to 
handle this overwhelmingly complex data to reveal the 
underlying physiological variations. Fourier analysis simply 
30 does not work for such nonstationary signals. By applying the 
invention, however, which works particularly well with 
nonstationary signals, one can arrive at a much deeper 
understanding of the underlying physiology. 
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As described above in detail, the invention applies a 
unique Sifting Process. Figures 15(e) -(h) illustrate this 
Sifting Process as it is applied to the blood pressure signal 
of Figure 15(b). Specifically, Figure 15(e) illustrates 
5 generating the envelopes (dot-dash line) by connecting 

successive extremas (peaks and troughs) of the signal (solid 
line) with cubic splines. The mean of these two envelopes is 
shown in Figure 15(e) with a thick solid line and is denoted 
by the symbol m 1 (t) . 

10 The difference, X(t)-m 1 (t) is designated h x (t) as is shown 

in Figure 15(f). From this figure, it can be seen that h x (t) 

□ has a few negative local maxima and positive local minima, 

suggesting that further Sifting needs to be performed before 

W the first intrinsic mode is generated by the invention and 

gf that the quantity h x (t) is not yet a true representation of an 
" s "oscillation about the mean." 

O To improve this situation, two definitive requirements 

ft j 

Hi are imposed for a function that represents the "oscillations 
W about the mean": (i) in the whole data set, the number of 

JSS* 

ftp extrema and the number of zero-crossings must either be equal 
or differ at most by one, and (ii) at any time, the mean value 
of the envelope of the local maxima and the envelope of he 
local minima must be zero. As oscillatory function that is 
processed to satisfy these definitions is called an intrinsic 

25 mode function (IMF) as more fully described above. 

The function h x (t) shown in Figure 15(f) does not satisfy 
this definition or the requirements of an IMF. Thus, further 
Sifting is performed by treating h x (t) as the new data set and 
repeating the Sifting Process including determining the new 

30 upper and lower envelopes, computing their new mean m xl (t) . 

The difference h 1 (t)-m 11 (t) is designated at h l:L (t) which is 
treated as the new data set. The process is repeated a number 
of times (Figs. 16 (g) - (h) ) until it converges. The convergent 
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result, Figure 16(g), is designated by C^t) and is the first 
intrinsic mode function, which has a zero local mean. 

The difference between X(t) and C x (t) is a function of 
time representing a "mean tread" after the first round of 
5 Sifting. It is designated as the "first residue" R^t): X(t)- 
C 1 (t)=R 1 (t) 

The first residue, R^t), is again oscillatory and can be 
analyzed as new data by another round of Sifting in the 
Empirical Mode Decomposition Process, yielding the second 
10 intrinsic mode function C 2 (t) and the second residue R 2 (t). 
M The process is continued until either the residue of the 
iZ\ intrinsic mode function becomes less than a predetermined 
small number of the residue becomes nonoscillatory . The 

UJ 

-dj resulting IMFs are shown in Figures 16 (a) - (p) . 

9.^ If the process takes n steps, the original signal X(t) 

s can be represented as follows: 



X( t) =C X ( t) +C 2 ( t) + . . . C n ( t) 24 



The last term C n (t) represents the nonoscillatory mean 
trend of the signal. As such, important information exists in 
20 this mean trend particularly for certain types of data. 

The other terms C^, C n _ 2 , etc represent the oscillatory 
portion of the mean trend. Each term represents a different 
spectral portion of the mean trend. 

With the IMFs in hand, further processing is possible. 
25 One such processing technique is filtering in which certain 
IMFs are combined or summed while leaving out other IMFs . 

For example, a low-pass filter or low-frequency 
representation M k (t) of the mean trend of the signal X(t) can 
be generated as follows: 

M k (t)=C k ^C k ^...C n 25 

30 
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where 2<k<n 

The lower the value of k, the more oscillations M k (t) 
contains. By adjusting the values of k and n, a variety of 
filters may be constructed each having desired 
5 characteristic (s) . 

Representing the original data, e.g. of Figure 15(a), 
with an analytic expression is not always possible or 
practicable. The invention, however, offers a technique of 
accurately representing the physically meaningful portions of 
10 a signal with an analytic function. 

For example, when one examines the indicial functions 
that represent the changes of the signal X(t) in response to 
step changes in oxygen concentration, one looks for changes in 

"3. S 

yj the mean treads M^t), M^Ct), and ... with respect to changes 
Jt£ in 0 2 level. To do this, it is helpful to fit the experimental 
S 1 result on M k (t) with an analytic function. Particularly, if 
p one takes the origin of time t=0 at the time of imposing a 

step change in oxygen concentration, the resulting change in 
the mean blood pressure in response thereto for t>0 can be 



a 



ill 



J?p represented as 



M k ( t) =A 1 +A 2 e" 1 2 fc +A 3 te"* 1 3 fc + A jn te" 1 mt 26 



where constants A 1# . . . and _ 2 , . . may be determined with 
25 a conventional least-squares estimation technique. 

The above estimation is greatly affected by the choices 
of k and m. By appropriately selecting k and m, different 
degrees of detail can be presented. 

30 In reality, one can only approximate a step function 

change of 0 2 . The dropping step of Figure 15(d) and the rising 
step of Figure 15(c) cause different changes in blood 
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pressure. Hence, two different empirical functions are 
utilized. For the first case, the equation is: 



fc-tp t-t 0 

— 1 T> - 7\ _ / i- \ 1" 



M k (t)=p m (t Q ) [l+A 1 (t-t 0 )e- ± T ± +A 2 (t-t 0 ) e" 1 T 2 

-l '-'o 

+A 3 [l-e 3 Ti ]], for t 0 < t (27) 
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for t 0 < t, where t 0 is the instant of time when 0 2 
concentration drops suddenly, and p m (t 0 ) is the measured value 
of M k (t 0 ). For the second case the following equation appears 
good enough : 

M /C (t)=p jn (t 1 ) [l+A 4 [e- 1(t - t l ) -l] ] , for t x < t, (29) 



where t 1 is the time when the 0 2 concentration increases 
suddenly. The variables A lf A 2 , A 3 and A 4 are dimensionless . 

Taking the analysis into the spectral domain, the 
invention offers further tools. As defined above, the Hilbert 
transform of X(t) is Y(t). Hilbert has shown that the complex 
variable Z(t) = X(t) + iY(t) is an analytic function of t and 
can be written in polar coordinates as a ( t ) exp { i_ ( t ) } , thus 
defining the amplitude a(t) and phase angle __(t). 

As further defined above, an instantaneous frequency _(t) 
can be calculated from the derivative of _(t) with respect to 
time. Knowing the amplitude a(t) and frequency _(t) as 
functions of time, one can plot contour maps as shown in 
Figure 23 (a) . The contour maps of the amplitude as functions 
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of frequency and time are called the Hilbert amplitude 
spectrum, H(_,t) . 

The vanishing of the local means of the IMF's C lf . . . C n ^ 
is a very important result because the amplitude a(t) and 
phase angle _ (t) of the Hilbert spectrum are sensitive to the 
local means . 

Using the definition above {M k (t)= C k + C k+1 + . . . C n , 
where 2<k<n} to represent the mean trend, the invention 
further defines the corresponding sum as follows: 

X A: (t)=C 1 (t)+C 2 (t)+. . .C k (t) 30 

to represent the oscillations about the mean trend M k+1 (t) . The 
Hilbert amplitude spectrum of X k (t) may be designated as 
H k (_,t). The square of H k (_,t) represents an oscillatory 
energy density. Thus, the invention also defines the 
oscillatory energy about the mean, E k (t), by an integration 
over all frequencies as follows: 

U 2 

E k (t)= w H k (w, t) dw 31 

The variations of E k (t) and M k (t) with oxygen level 0 2 (t) yield 
the desired summary of information about the transient 
response of X(t) to 0 2 (t). 

The oscillatory energy about the mean E k (t) is quite 
similar to the instantaneous energy instantaneous energy 
density, IE, described above. The difference is that the E k (t) 
is calculated from a subset (k) of the IMFs while IE is 
calculated from the full set of IMFs. In other words, the 
oscillatory energy about the mean is derived from a filtered 
version of the original signal in which a subset (k) of the 
IMFs is utilized to construct the filtered signal. This 
contrasts with instantaneous energy density which is derived 
from the full set of IMFs 
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The oscillatory energy about the mean E k (t) is another 
nonstationary stochastic variable that can be treated in the 
same manner as outlined above. Specifically, the progression 
from steps 350 to 370 to 120 in Figure 1(d) provide a process 
5 for applying EMD to E k (t) to generate another set of IMFs . The 
IMFs thus generated may be further processed according to the 
techniques disclosed herein. The process may be repeated as 
many times as desired as indicated by the loop in Figure 1(d) . 

From the IMF's shown in Figures 16 (a) - (p) , one can 
10 compute the mean trend functions M k from equation 25 set forth 
above. The results are shown in Figures 17 (a) -(f) which are 

ll[ the ma -i n results to be used for the indicial response 

i J 

Q determination . 

-••A 

u\ The next step (360 Figure 1(d)) is to generate analytic 

Cfe functions or otherwise perform a curve fitting process. This 

in 

%] is a chief advantage of the present invention because for many 
J data sets, such as the blood pressure data in Figs. 15 (a) -(d), 
deriving analytic functions is impossible. By applying the 
inventive methodologies, however, deriving such analytic 
pb functions is possible. Specifically, the mean trend data 

extracts relevant data and removes unnecessary or irrelevant 
data thereby producing a simplified data set for which it is 
possible to derive an analytic function. The degree of 
simplification and extraction can be controlled by properly 
25 selecting which IMF components are included in the mean trend 
data . 

The analytic representation of the indicial response 
functions of Figures 17 (a) -(f) by equations 26-28 was done 
with a least-squares error method. The results are shown in 
30 Figures 18 (a) -(c) and 19 (a) -(c). These are the primary 
results of interest for tissue engineering analysis and 
design . 

Specifically, the results shown in Figures 18 (a) -(c) and 
19 (a) -(c) provide indicial response functions of the mean 
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pulmonary blood pressure at the arterial trunk in response to 
a step decrease or a step increase in breathing gas oxygen 
concentration. A similar analysis can be directed toward the 
oscillation modes defined by equation 29. 
5 The oscillation modes for k=l, 2, 4, and 6 are 

respectively shown in Figures 20 (a) -(d). The signals shown in 
these figures contain a lot of information that needs to be 
extracted into simple, understandable terms. This can be done 
by applying the Hilbert transform and generating the Hilbert 
10 spectrum. The results are shown in Figures 21, 22 (a) -(h), 
H 23 (a) -(b). Specifically, Figure 21 is a plot of the 
p oscillatory energy defined by equation 29 as a function of 
time. This spectrum is analyzed as a nonstationary random 
v=j signal by the invention to resolve it into oscillatory IMF 

m 

-If modes and mean trend functions M k (t) , with the associated 

= analytic functional representation as illustrated in Figures 

fjj 22 (a) -(h) for k=9, 10, and 16. 

\'\ To further analyze the data, the results of calculating 

p the instantaneous frequency and amplitude of oscillations of 
s 2t> the pressure signal X(t), made possible by the Hilbert 
transformation, can be plotted against time in a three- 
dimensional manner as sown in Figure 23 (a) . The same can be 
plotted two-dimensionally on the plane of time and frequency 
by the contour lines of equal amplitudes, as show in Figure 
25 23 (b) . 

Figure 24 is a conventional FFT amplitude spectrum of the 
pressure signal in 1-min segments under the assumption that 
the process is stationary in each segment. This conventional 
spectrum provides a basis of comparison against the Hilbert 
30 spectra generated by the invention. 

As mentioned above, Fourier analysis is based on the 
hypothesis of segmental stationary random oscillation, the 
principle of linear superposition of sine waves, and a global 
average of waveform convolution over each time segment. The 
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HHT is based on the hypothesis of nonstationary processes, the 
principle of linear superposition of nonlinear IMF's, and 
local determinations of amplitude and frequency (through 
differentiation rather than convolution) of each IMF. In 
terms of the IMF, the first k modes can be added together the 
represent oscillations about the mean trend M k+1 (t) . 

The Fourier series cannot represent time variation of a 
nonstationary signal and does not have a property or 
capability to separate a signal into two parts, one part 
representing a mean trend while the other part represents 
oscillations about the mean. The number of intrinsic modes, 
n, is finite: in general, n<log 2 N, where N is the total number 
of data points . The number of harmonics in Fourier analysis 
is N/2 . 

Comparison of the Hilbert and Fourier spectra 
respectively shown in Figures 23 (b) and 24 shows that 
both spectra display a major frequency component at 
approximate 5Hz where the energy is concentrated. This is 
close to the heart rate of the rat subject. This rate 
decreases when the oxygen concentration is decreased. These 
two spectra are in different vertical scales. The Hilbert 
spectrum contains no energy with frequency > 10Hz, and it also 
has fewer yet more diffused frequency bands that the Fourier 
spectrum. This is because the Hilbert spectrum gives only the 
global mean. The mean values certainly will show less 
variations . 

The Fourier spectrum of Figure 2 4 contains more frequency 
bands because any deviation of the waveform from the basic 
sinusoidal harmonic will result in strong higher harmonics 
whereas the Hilbert spectrum allows variation of instantaneous 
frequencies, hence the fuzzy spread in the frequencies. This 
calls attention to the fact that the heart rate is also a 
stochastic variable, which could be studied by the inventive 
methods disclosed herein. The strong higher harmonic band 
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with frequency > 10Hz in the Fourier spectrum is probably 
spurious. In other words, the Fourier components above 10Hz 
(staring at around 15Hz at time 0) have no physical meaning. 
In contrast, the inventive results of Fig. 23(b) show a 
5 rich variety of low frequency components thereby conveying 

important information not available in the conventional result 
of Fig. 24. Thus, the invention is capable of more a more 
accurate physical representation of the underlying phenomenon 
than conventional Fourier-based analysis. Moreover, the 
10 clarity of the set of mean trends and the corresponding set of 

the oscillations is a unique contribution of the invention. 
III Thus, the invention provides useful tools for concisely 

O and precisely describing the affect of changing one variable 
!„= on another variable. The effect on blood pressure caused by 
€k changing the 0 2 concentration provides a good illustrative 

in 

example of these tools . 

To further illustrate the broad applicability of the 

Q 

fU invention and its potential as a diagnostic tool, an abnormal 

L=i 

: , condition (disease) was chosen. The particular abnormal 

*2D condition chosen for illustration purposes is sleep apnea. 

This is a common condition in which the airway is temporarily 
obstructed during sleep causing the subject to eventually gasp 
for air. To study this condition, heart rate data was 

taken from a subject as shown in Figure 25. More accurately, 

25 the data is a measure of the time interval between pulses 
(pulse interval) of the heart plotted against time. Thus, 
small values indicate a fast heart rate and large values 
indicate a slow heart rate. This data includes a normal data 
section and an abnormal data section, both of which are 

30 labelled. 

As will be shown in detail below, the invention provides 
powerful tools for studying this condition. Although the 
results to date are preliminary, they do indicate that the 
invention is capable of diagnosing this condition. Of course, 
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further studies are necessary before a reliable diagnosis can 
be made using the invention. 

To improve the analysis, the data of Figure 25 was 
processed to make the spacing between data points even. 
5 Furthermore, a section of the normal data (no apnea episode) 
was extracted for separate analysis. A blowup of the normal 
data having even spacing between data points is shown in 
Figure 26. A further blowup is shown in Figure 28. 

The data of Figure 2 6 was then subjected to the inventive 
10 EMD to produce eight IMFs, C^Cq, as shown in Figures 27(a)- 

(h) . The resulting IMFs provide important, and heretofore 
r .x unavailable, data to a person attempting to gain a better 

understanding of sleep apnea particularly when these normal 
y IMFs are compared against abnormal (apnea) IMFs. 

S3. 

'is The IMFs of Figures 27 (a) -(h) were then utilized to 

in 

Sj construct the Hilbert Spectrum shown in Figure 29. This 
m Hilbert Spectrum provides additional information about normal 
fll sleep patterns and how they are reflected in the subject's 
Lj heart rate. 

jito To continue the analysis, an abnormal data section 

(Figure 30) was subjected to the inventive EMD to produce 
eight IMFs, C^-Cg, as shown in Figures 31(a)-(h). The third 
IMF C 3 is believed to be particularly important in analyzing 
and perhaps even diagnosing sleep apnea. 

25 Figure 32 is a blowup of the abnormal data of Figure 30. 

The abnormality is apparent around second 5005 where the 
airway is temporarily blocked thereby causing the heart to 
race in an effort to extract the diminishing oxygen still 
available in the lungs. The sharp spike just after second 

30 5005 illustrates the gasping breath and the resulting sharp 
increase in the pulse rate interval (corresponding to a 
decrease in the heart rate) . 

Figure 33 is the resulting Hilbert Spectrum of the 
abnormal data of Figure 30. Comparisons between the normal 
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Hilbert Spectrum of Figure 3 9 with the abnormal Hilbert 
Spectrum of Figure 3 0 yields further information for studying 
and understanding sleep apnea. 

To further illustrate the wide applicability of the 
5 invention, it was applied to study epilepsy. Epileptic 

seizures occur when the brain's neurons fire in synchronism. 
The sequence starts from an epicenter and then propagates to 
an entire hemisphere of the brain. Even at this stage, the 
patient can still function in a relatively normal manner. 

10 When the synchronization propagates to both hemispheres of the 
brain, then the patient will suffer a seizure. 

f=j To study epilepsy, heart beat rate data was chosen. 

t| Figure 3 4 shows heart rate data (beats per minute) that was 

U\ collected from a patient before, during and after a seizure. 

£5 The time 0 is the starting point of the seizure. 

%l Figures 35 (a) -(i) show the IMF generated by the inventive 

~! EMD method applied to the Figure 34 data. These IMFs reveal 
fU important and physiologically meaningful data not available to 

sss's 

hi a researcher before the application of the present invention. 

20 Figures 36(a) and (c) illustrate two IMFs (modes 1 and 3) 

ill 

which may be utilized to correlate the data and study 
epilepsy. These IMFs are shown in alignment with the original 
heart rate data (Figure 36(b)). As is apparent from these 
figures, the first and third have components that appear to 
25 correlate with the onset of a seizure. 

The Hilbert Spectrum of Figure 37 was generated by the 
invention from the data of Figure 34. The trace along 18 to 
20 cycles/min correlates to time 0, the start of the seizure. 
To compare the invention against conventional techniques, a 
30 corresponding Wavelet Spectrum is shown in Figure 38. 

As can be seen by the above illustrative examples, the 
invention is a powerful tool capable of analyzing data and 
producing new, quantitative measurements that enhance the 
understanding of the underlying phenomenon producing the data. 
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The IMFs themselves are important results of the inventive 
data analysis techniques and permit one to gain a level of 
understanding not possible with conventional techniques . 

Although blood pressure, heart pulse interval, and heart 
5 rate data provide good examples of the invention's 

applicability to biological signals, it is to be understood 
that a large variety of other biological signals may be 
processed by the invention to gain a better understanding of 
the underlying biological process (es). Other examples include 
10 plethysmogram signals, electro-encephalogram (EEG) signals and 

temperature signals. Furthermore, the signal need not be 
ll[ taken from a living body and include in vitro studies such as 
ill current flow across membranes (patch clamping) , fluorescence 
|.l in confocal microscopy, spectroscopic signals, etc. 

tfs As mentioned above, the invention is also not limited to 

\ n 

biological signal processing and includes the full range of 

1^ real-world, data representative of an electrical, chemical, 

O 

f|j mechanical, optical, geophysical processes or phenomenon or 

n combinations thereof all of which fall under the term 

Kb "physical signal" as it was defined in the parent application. 

Even further, the invention provides tools for studying 
the influence of one variable on another variable in a multi- 
variable system. The effect of hypoxia on blood pressure 
provides one illustrative example of this analysis. However, 
25 it is to be understood that any physical signal in a process 
or phenomenon having multiple variables of which the physical 
signal is one can be analyzed with the invention. 

ALTERNATIVE EMBODIMENTS OF BIOLOGICAL SIGNAL ANALYSIS 

30 As described above, the invention constructs upper and 

lower envelopes 20, 30 with a cubic spline in steps 210 and 
230, respectively and in step 560. This cubic spline fitting, 
however, has both overshoot and undershoot problems. These 
problems can be alleviated by using sore sophisticated spline 
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methods, such as the taut spline in which the tension of the 
spline curve can be adjusted. 

Another alternative is higher-order spline fitting. 
Although such higher-order spline fitting may be more 
5 accurate, it will, however, introduce more inflection points 
or extrema, and consume more computing time. Therefore, it is 
not recommended as a standard operation. Only in special 
cases, it may be tried. 

As the spline fitting procedure is time consuming, more 
10 efficient methods can be devised by using simple mean values 
M of successive extrema instead of computing the envelope-mean. 
pi In this way, only one spline fitting is required rather than 

N two. Although this alternative is easier and faster to 

W 

,Pj implement, the shortcomings are more severe amplitude 

%h averaging effects when the neighboring extrema are of 

s different magnitudes. The successive-mean method will have a 

?~% 

?H stronger forcing to reach uniform amplitudes, in which the 

i ~-j 

H true physics associated with amplitude will be destroyed. 

Ul 

f =l Therefore, the successive-mean method should only be applied 
6b where the amplitudes of the physical signal components are 
constants . 

Either the envelope mean or the successive mean method, 
when applied with the requirement of absolute symmetry, will 
produce the absurd result of uniform amplitude IMF's. 

25 Therefore, the criteria in the Sifting Process should be 
chosen judiciously. One should avoid too stringent a 
criterion that we would get uniform amplitude IMF's. On the 
other hand, one should also avoid too loose a criterion that 
would produce components deviating too much from IMF's. 

30 It is well known that the most serious problem of spline 

fitting is at the ends, where cubic splines can have wide 
swings if left unattended. As an alternative, the invention 
may utilize a method of adding characteristic waves at the 
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ends of the data span. This confines the large swings 
successfully . 

The method of adding characteristic waves is not 
conventional. In contrast, the conventional window that is 
often applied to Fourier transform data results in loss of 
useful data. To avoid this data loss and to confine swings at 
the ends of the data span, the invention extends the data 
beyond the actual data span by adding three additional 
characteristic waves at each end of the data span. 

The characteristic waves are defined by the last wave 
within the data span at the end of the data span. In other 
words, a characteristic wave is added to each end of the data 
span having an amplitude and period matching the last wave 
within the data span. This characteristic wave is a 
sinusoidal waveform that is extended three sinusoidal wave 
periods beyond the data span at each end. This process is 
repeated at the other end of the data span. In this way, 
spline fitting at the end of the data span, which can 
otherwise have a wide swing, is confined. In other words, by 
adding the extra characteristic waves at the ends beyond the 
data span, the spline curve will be tied down so that it will 
not have wild or excessive swings that would otherwise corrupt 
the data processing and analysis that utilizes these cubic 
splines . 

Other than the spline fitting, the Hilbert transform may 
also have end effects. Because the first and the last points 
of the data are usually of different values, the Fourier 
transform will introduce additional components to bridge over 
the difference resulting in the well-known Gibbs phenomena. 
To eliminate it in the Fourier transform, various windows have 
been adopted (see, for example, Brigham, 1974, "The fast 
Fourier Transform", Prentice-Hall, Englewood Cliffs, NJ) . 

Instead of a window which will eliminate some useful data 
at the end, the invention again adds two characteristic waves 
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at either end. These waves all start from zero at the 
beginning, and end at zero at the end. Thus, the annoying 
Gibbs phenomena are greatly reduced. 

Still further, the Hilbert transform needs over-sampled 
5 data to define the instantaneous frequency precisely. In 
Fourier analysis, the Nyquist frequency is defined by two 
points per wave, but the frequency is defined for a wave 
covering the whole data span. In the invention, the 
instantaneous frequency is defined through a differentiation 
10 process, and thus more data points will help defining the 
frequency more accurately. Based on the inventor's 

U experience, a minimum number of data points to define a 

Q 

G frequency is five (or four _t ' s ). The lack of fine time 
4; steps can be alleviated through interpolating more points 
m> between available data. As a spline interpretation would also 
not create nor annihilate scales, it can also be used for the 

(3 interpolation when the data is very jagged from under- sampled 

rl 1 

[2 data. The smoothed data though have a longer length and are 

M sometimes easier to process. The interpolation may give 

Hp better frequency definition. 

ACOUSTICAL SIGNAL ANALYSIS 

Referring to Figure 39, an application of the invention 
will be to analyze the acoustical signal through the Hilbert- 

25 Huang Transformation (HHT) , which consists of the Empirical 
Decomposition Method (EMD) and the Hilbert Spectral Analysis 
(HSA) . Essentially, the signal will be decomposed into the 
Intrinsic Mode Function Components (IMFs) . This decomposition 
is outlined in Huang's publications and patents. Once the 

30 invention decomposes the acoustic signal into its constituting 
components, all the operations can be performed on these 
components . 

Speech analysis and speaker identification: 
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Referring to Figures 41 and 42, we will first discuss the 
way HHT analyzes speech. Then, based on the analysis results, 
we can discuss how to identify a speaker, i. e., to single out 
one speaker from others. Let us start with the following 
5 example. 

Two speakers both said 'Halloo, ' and the acoustic signals 
are digitalized at 44100 Hz with 16 bits resolution as shown 
in Figure 45a for a male speaker, and b for a female speaker. 
The EMD decomposed IMF components are given in Figure 4 6a and 
10 b. Each signal consists of 12 components, and each component 

consists of 30,000 points. They are plotted with the same 
U vertical scale in Figure 46. As the IMFs are the results of 
j==^ scale separation of the original signals, each component 
y represents different time scales, which can be translated to 

^5 the frequency ranges. It should be pointed out that this 

y i 

%| frequency separation is not in the Fourier sense, for the 

f-% signals separated in scale could be nonlinear; therefore, it 

III can contain both sub and super harmonic components in Fourier 

u\ frequency space. Consequently, the decomposition is not 

Vp necessarily narrow banded in Fourier sense, but, as explained 

fl? 

in Huang (1999), this decomposition preserves all the 
nonlinear properties of the signal. One can immediately see 
from these figures that the two decompositions contain very 
different IMFs. For the signal given in Figure 45a, the 

25 decomposition is given in Figure 46a. There, we have very 

regularly periodic signals for all the IMFs. Comparing this 
with the corresponding results in Figure 46b from the second 
speaker, we can see the differences. The second speaker does 
not have the regularity in the time scale. In the second 

30 speaker's results, we can see that the strongest signal is the 
third component, but for the first speaker the corresponding 
component is relatively weak. By the time we come to the 6 th 
component of the first speaker, the component is still very 
strong, but there is almost no more signal in the second 
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speaker's result in Figure 46b. The richness of the low 
frequency components for the male speaker, and high frequency 
for the female speaker should be expected. This difference in 
frequency components can be seen vividly in the Hilbert 
5 spectral analysis as follows. 

From these components, we can construct the Hilbert 
Spectra as given in Figures 47a and b, in which time- 
frequency-energy representation of the signals are presented, 
with the corresponding narrow banded Fourier based 
10 spectrograms given in Figures 48a and b constructed with 102 8 
Ll sampling points in the window, which corresponding to the 

r 

O window width of 43Hz. In these figure, the total frequency 

f "i 

range is represented by 2 00 bins to cover the frequency range 
'I? from 0 to 10,000 Hz. In Figure 47a, we can identify a clear 
g5 consistent low frequency energy concentration near 100 to 200 

Hz, but there is no equivalent energy at this low a frequency 

U in Figure 47b. To examine the results in more details, we 

ry 

12 have given the Hilbert spectra of the cases presented in 

Figure 47a and b in Figures 49a and b, in which the same 2 00 

np frequency bins are reassigned to cover the frequency range 
from 0 to 1,000 Hz. In these new and more detailed 
presentations of the signals, we can see that the male voice 
represented by the Hilbert spectrum by Figure 49a contains a 
strong modulated energy below or around 100 Hz, while the 

25 female voice represented by the Hilbert Spectrum given in 

Figure 49b contains the most energetic modulated components 
above 100 Hz. Furthermore, the female voice also contains 
wider range of high frequency energy distributions. In Figure 
49, we can see the regular spaced vertical striations in the 

30 Hilbert spectra, which are the indications of the resonant 
beat period from the speaker's vocal cavities. This beat 
pattern can be easily seen from the IMF components, which can 
serve as crucial data for speaker identification too. 
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In order to show the unique properties of the HHT 
analysis results, we decided to examine the male voice in more 
details as follows: First, we give the spectrum it wide 
banded Fourier spectrogram as in Figure 50, constructed from 
5 12 8 sampled points in a window, which corresponds to a wide 
banded spectrogram with band width of 344.5 Hz. Clearly 
neither the narrow nor the wide banded spectrograms give us 
any detailed information. This is shown even more clearly if 
we present the detailed views in Figure 51a and b, in which 
10 the spectrograms are plotted for the frequency range between 0 
M to 1000 Hz, corresponding to the condition represented in 
yd Figure 49a. The comparison is clear, while the Fourier based 

%I spectrograms could give us some qualitative information, the 

til 

,fj results is so crude that they cannot be used as a based for 

•-ife any serious acoustic signal analysis. 

s In order to emphasize the superior performance of the 

JjJ HHT, we also process the same data with Morlet Wavelet 

h= analysis, a continuous Wavelet analysis. The result is given 

W 

?* % in Figure 52a and b. These figures clearly identify the 

lib different time-energy- frequency distribution characteristics 
between male and female speakers, the strong high frequency 
energy concentration for female and the strong low frequency 
energy concentration for the male speaker. This clear 
difference aside, the Wavelet results really lacks the 

25 detailed analysis ability of HHT. When on comparing the pair 
of Figures 52a and b with those in Figures 49a and b, one can 
immediately see the lack of time- frequency resolution of the 
Wavelet results. These differences have been discussed in 
great details by Huang et al (1998, 1999) . 

30 To illustrate the lack of resolution in the Wavelet 

analysis further, we can smooth the Hilbert Spectrum result 
given in Figure 49a by filters in 100 and 200 time-steps. The 
results are given in Figures 9a and b respectively. By the 
time we smoothed the results with filter width of 200 points, 
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the results are clearly much smoother, but we will never 
obtained the kind of overly smoothed result given by the 
Wavelet, for the two analysis methods are very different: The 
Hilbert spectrum gives the genuine instantaneous frequency 
based on an adaptive basis totally independent of the kind of 
a priori basis adopted in the Wavelet analysis. Furthermore, 
the components in HHT are almost orthogonal, which gives no 
energy leakage among different modes. The continuous Wavelet, 
on the other hand, has a non-orthogonal representation. The 
energy leakage among different frequency bands, and the over 
redundancy have rendered the method to only extract 
qualitative information . 

The detailed time-frequency-energy presentation gives us 
many more reference points for speaker identification and 
verification. The continuous modulation of frequency and 
energy will provide more reference points to separate one 
speaker from the other. HHT analysis could also identify 
peculiarities of one speaker from another. All these details 
just cannot be found in any other analysis method. 

To summarize the material discussed in this section, we 
will first give a flow chat of speech analysis: Speech 
containing various phonemes will be recorded, and the signal 
will pass through the sifting process to get the IMFs . In 
this step, the critical task is to identify the phonemes and 
formants (Hilbert Spectra) for each sound element with 
modulations in both frequency and amplitude as functions of 
time. Databases will be build from the recorded speech data. 
Such analysis is not available from any other method of 
analysis . 

There are two types of speaker identifications: to 
verify that the speaker is a known person, and to separate the 
speech from different speakers. The essential steps for both 
tasks are as follows: 
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To verify a known speaker, we have to have the database 
of the targeted speaker built ahead of time. The database can 
be built following the steps listed in the speech analysis. 
To separate speakers from one another, we can record the 
5 speakers and analyze their speech as listed in the speech 

analysis part. Once the analyses are done, we can compare the 
details of the Hilbert spectra in their over all and the 
details of the spectral patterns, including the beat periods 
for the same phoneme and formant. Given the super temporal 

10 and frequency resolutions of the Hilbert spectra, the Hilbert 
spectra actually give a much better formant representation. 
With so many details in the Hilbert spectra, we can use many 

p features as references. In the example given in this section, 
the overall pattern is sufficient. 

in 

: ; For Speech recognition: 

= Speech analysis such as speech recognition as shown in 

r|j Figure 42, the goal is not to identify a specific speaker, 
f" Rather the task is to identify the meaning of an utterance so 
Sp that the commands contained in the speech can be translated 

s: : 
: s ; 

5 " into actions as in man-machine interaction. Because, the 
speakers can have a variety of accents, the approach is to 
obtain a limited number of uniquely identifying parameters to 
define the meaning of the utterance, rather than to obtain 

25 more and more detailed descriptions of the sound. In such 
applications, two approaches can be adopted for the HHT . 

The first is by smoothing as in the cases of Figure 53a 
and b. More smoothing can be obtained by time- frequency 
averaging, rather than just time smoothing described here. 

30 The Laplacian filter is readily available. However, for most 
acoustic signal applications, there are lot more time steps 
than possible frequency bins. Therefore, time-wise smoothing 
is a more logical approach. Other than smoothing to obtain 
the time- frequency-energy pattern as means for recognition, 
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the formants identification is also a powerful method. Here 
we will discuss this approach as follows: 

Formants by definition are the few characteristic time- 
frequency distribution of the maximum energy concentration in 
a speech section analyzed with filter bank or spectrogram 
methods (see, for example, Rabiner and Schafer, 197 8) . 
Conventionally, to make the number of components manageable, 
the formants are defined from the low frequency components 
only. The ideal of formant is to sharpen the poor time- 
frequency analysis in the spectrograms. Similar approaches 
have also been explored by Carmona et al (1997) using Wavelet 
analysis. As discussed by Huang et al (1998, 1999), the 
inherited poor time- frequency resolution cannot be overcome by 
post analysis smoothing or selection. The only alternative is 
to change to a genuine instantaneous frequency analysis as in 
the HHT, which in form and in substance is a much better 
formant representation of the speech. Let us examine an 
example. If we take the beginning of the sound for 'Halloo' 
signal, and analyze it with HHT, the result is given in 
Figures 54a and b. This is actually the consonant part of the 
x h' sound. These figures, in fact, are simply obtained by the 
enlargement of the same section of the figure of 47a and 49a. 
One of the great advantages of HHT is that the frequency is 
defined by differentiation; therefore, we can breakthrough the 
limitation of the uncertainty principle and expanding the 
results to the limit of the sampling rate. 

Figure 54a is the expansion of Figure 47a, while Figure 
54b is the expansion of Figure 49a, both with the original 
signal plotted. Here we have the lower frequency emphasized. 
The continuous lines of the time- frequency curve are in form 
and substance a better formant. But they are indeed better 
defined than the conventional formants, for they also included 
the nonlinear deformation of the signal rather than assigning 
all the waveform deformation to harmonics. As the frequency 
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is defined by differentiation, the information provided here 
is might be too much. If that is the case, the frequency- 
variation can be simplified by smoothing either through 
filtering, or defined as differentiation with larger time 
5 steps. 

The essential steps for this task can be summarized in 
the following flow chart: Again, we have to establish the 
database consisting of templates of the key phonemes and 
formants. In this case, the templates will be constructed 
10 from the averaged value of a large number of different 
Ll speakers. In addition to the averaged value, the range of 

O scattering from the differences will also be recorded. When a 

C3 

\l new speech is given, the record will be analyzed as in the 
\!; previous part to have the speech breakdown to phonemes and 
tfb formants. Each sound element will be compared with the 

existing templates and identified. The identified sound will 

G be put together, and translated into sentences. Grammar check 

ft! 

l\ and corrections should be made here to guarantee that the 

speech makes sense. The final text will then be translated 

SO into instruction for further actions, such as an order to a 
servomechani sm . 

Speech synthesis: 

Referring to Figure 40 the speech synthesis problem has 
25 been review recently by Breen (1992). Basically, the 

synthesis is based on a phoneme base. Having established the 
analysis part, we can take sections of the sound record as 
phonemes and build a database for it. This phoneme database 
can then be used to synthesize speech. The advantage of this 
30 phoneme base has the unique advantage over the conventional 
one in that no local stationary assumption has been evoked. 
Therefore, the frequency is continuously modulating, which can 
make the final result sounds more naturally. 
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The detailed time- frequency-energy distribution results 
given by the HHT also give us a means to synthesize sound with 
higher fidelity, for the frequency modulation can be truly 
simulated by the instantaneous values in frequency and energy. 
5 This overcomes one of the fundamental shortcomings of the 

Fourier based analysis, in which a local stationary assumption 
is to be invoked. This piecewise stationary analysis simply 
cannot closely simulate the constantly modulated frequency 
modulation as shown in the Hilbert spectral analysis results. 
10 Also from Figures 54a and b, we can define the beginning 

L = 3 of the speech clearly and un- ambiguous ly . The separation of 
J;J syllables has also been a hard problem in speech analysis and 
S| synthesis. 

a _ ; 

l}\ A key to the potential value of our approach is that HHT 

tfl gives time varying modulation of frequency and amplitude. The 

s essential steps are as follows: 

O 

32 5 
111 

M Based on the database of all the phonemes and formants, we can 

til 

translate the given text to the time varying IMF components 
Ub from the database. This translated data is checked again with 
prosody patterns, and synthesized into sound. 

Speech quality enhancement: 

Referring to Figure 43 one of the most important 
25 applications of acoustic signal analysis is to enhance the 

speech quality. For example, during the utterance of the word, 
x Halloo, ' a noise simulated by the drop of a spoon on the 
table was also heard. How to separate this noise from the 
speech signal is an important application. Traditionally, the 
30 cleaning up of a speech is carried through filtering But 

Fourier based filtering has inherited shortcomings: Localized 
signal in time is broad in frequency; localized signal in 
frequency is broad in time. Furthermore, as discussed above, 
the generating mechanisms of most consonant sound are 
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nonlinear; therefore, most of the consonants, the information 
carrying sounds, have a lot of sub as well as super harmonics. 
We cannot remove certain frequency without altering the 
harmonics of other sound signal, which will down-grade the 
5 filtered result through loss of the harmonics. Here we 

propose the filtering through EMD. The filtering strictly by 
scale, the filtering here has to be implemented jointly in 
time- frequency space . 

Let us examining the signal of the sound, x Halloo + 

10 ding,' given as the top line in Figure 55. This signal is 

decomposed through EMD and we obtained the IMF components as 

f=% given in Figure 5 6a and b. Figure 5 6a gives all the 

£3 components with the uniform vertical scale, while Figure 5 6b 

i 

~= 

U| gives the same information with the vertical scale in each IMF 

Si 

£5 components normalized within its own frame. From Figure 56a, 

in 

xl we can see that the x ding' sound, occurring slight after 1.5 
r*\ second mark, has the higher amplitude than any other IMF 

components. Figure 56b indicates that at the same time scale 
lil components with the x ding' sound, there are also signals from 
2§ the 'Halloo, ' as indicated by the low level signal in c3 . As 
a result, we cannot simply filter out the high frequency 
components by eliminating the first three IMFs . 

Now we will accesses the performance of the EMD vs 
Fourier filtering. First, let us implement only the filtering 
25 by simply eliminate the first three IMF components. This 

approach is similar to the thinking of Fourier filtering as 
discussed by Huang (2000) , but the result is not exactly the 
same. We will design this filtering as the "Fourier" type and 
discuss the difference later. The result certainly does not 
30 show the x ding' sound any more. But this brute force 

filtering has also eliminated all the high frequency sound 
including some higher frequency elements from the x Halloo' 
sound. The quality of the sound is degraded considerably. A 
better way is to eliminate ' the 'ding' sound in both time and 
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frequency space by selectively removing energy associated with 
the 'ding' sound. This was achieved by first studying each of 
the IMF component, identifying the 'ding' sections in each 
component by their sounds, frequencies, and occurrence time 
5 location. Then, selectively removing the sections that 

contain the 'ding' associated signals from each IMF component. 
After this is done, and separated 'Halloo' and 'ding' signals 
are represented by their respective IMF decompositions in 
Figures 57a and b. In Figure 57a, we can see that there is no 
10 longer signal associated with the 'ding' sound, while, in 
M Figure 57b, we see the clear 'ding' sound stands by itself. 
n The Hilbert spectra for the original and the EMD filtered 
'/\ signals are shown in Figure 58a and b. When all the IMF's are 
% i) summed up, the re-constituted sound signals for 'Halloo' and 
i ? 5 'ding' are given in Figure 55. The sound quality is so much 
s improved over the brute force approach. 

?H Next, we will implement the Fourier based filtering with 

H two 48 degree FIR filters with cutoff frequency set at 2 000 

h ! 

p and 1000 Hz respectively. The filtered data with a cutoff 

^6 frequency at 2 0 00 Hz together with the HHT filtered data is 
given in Figure 59a. Based on the results plotted here, one 
can hardly see the difference. But if one pays attention to 
the relative darkness of the two side-by-side signals, one 
would see that Hilbert filtered signal is darker, an 

2 5 indication of more high frequency components being retained. 
Furthermore, if one plays the filtered data back through a 
speaker, we can immediately tell the difference, for there are 
still residual 'ding' sound in the Fourier filtered case. We 
then filtered the data with a similar filter but with a cutoff 

30 frequency set at 1000 Hz. This certainly eliminated the 
'ding' sound, but the quality of the sound was very much 
downgraded . 

The Fourier filtering is strictly in frequency space, 
while EMD filtering is in time scale. By operating in a rigid 
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frequency range, the Fourier filter will not eliminate the 
bounded super and sub harmonics of the signal outside of the 
pre-selected filtering frequency range. It, however, will 
remove the bounded super and sub harmonics of other signal 
5 components that happened to fall within this pre-selected 
frequency range, even if the fundamentals are outside this 
frequency range. As a result, through Fourier filter is 
localized in frequency, it actually affects signals of all 
frequency. In EMD, each IMF could be nonlinear; therefore, 

10 each of them is not banded limited in Fourier sense. By 

eliminating the first three IMF components, we still could 

p have energy at Nyquist frequency in Fourier sense through the 

; s ! existence of the harmonics. So, through filtering by EMD, we 

""-4 

UJ can really limit the frequency to the components in time scale 

sa 

gjfjp rather then m Fourier sense. This time space filtering will 
not influence the lower frequency components because it would 

5 

p leave the bounded super and sub harmonics of other components 

r{i 

:* r untouched, and it will also eliminate all the bounded super 

?« 

IU and sub harmonics associated with the selected time scale 

ill 

Wo component. We have plotted the results of the Fourier spectra 

i u 

of the EMDT filtered, the "Fourier" type of filtering by 
eliminating the first three IMF components, the Fourier based 
FIR with cutoff frequency set at 2000 and 1000 Hz all in 
Figure 59b. As we discussed above, only by setting the cutoff 

25 frequency at 1000 Hz can we eliminate the x ding' sound. The 
spectra showed the difference: Fourier filtering has 
eliminated all the energy at the high frequency range, but EMD 
filtering can leave some energy there. The result is a 
greatly improved voice quality. 

30 Another application for the present invention is on music 

recording. Using EMD rather then Fourier based filter we can 
surgically remove the unwanted noises or sound to enhance the 
final quality of the product. As discussed above, the EMD 
filter is much more sophisticated than the Fourier based one, 
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for the EMD filter operates with the signal separated bi- 
temporal scales; therefore, it can keep track of the events 
precisely on the temporal axis . This operational principle is 
totally different from the Fourier based filter, which 
5 operates in Frequency space only. As discussed above, the 
localized signal in frequency space means the signal is 
uniformly distributed in time space, and localized signal in 
time space means uniformly distribution in frequency space, it 
is just impossible to perform the precision surgical-like 
10 filter localized both in time and frequency spaces. 

Furthermore, many of the sound full of harmonics will be 
O degraded no matter which frequency the Fourier filter will 
rt remove. In the EMD filter, however, the IMF components 
id retained all the harmonics within a single scale component, 
fife The harmonics are represented by the frequency modulation as 
% '4 discussed by Huang et al (1998, 1999) and Huang (1999) . 
p Consequently, any filtering operation will not downgrade 
IV signal at other scales, which can still retain all its 

rsss 

LiJ harmonics, if viewed in theFourier sense. 

3p A step in this task is to identify the unwanted sound in 

its time coordinate. It is impossible to blindly apply this 
method as in the traditional Fourier analysis to remove high 
frequency components. If that is the case, then we can also 
use EMD to remove the highest frequency IMF component. This 

25 approach has one unique advantage: By removing the highest 

frequency component locally or globally, it will not alter the 
quality of the low frequency sounds, for HHT does not depend 
on harmonics to represent nonlinear distortion of the 
waveforms of sounds. With the additional time tag, we can 

30 also implement HHT filter locally and globally with the 
precision unmatched by Fourier window-based filters. 

If the unwanted sound is identified on the time axis, we 
can zero in the IMFs and test each component, or a combination 
of several components, to determine which one or which 
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components should be removed. Once the unwanted component (or 
components) is removed, the cleansed sound could be re- 
constituted locally . 

Acoustical signal storage and transmission: 

Acoustical signal storage and transmission might require 
limited bandwidth or data quantity. Using HHT, one could also 
reduce the amount of data to represent a given signal with 
acceptable, albeit reduced, quality. This can be achieved by 
re-sampling of the IMF components of longer time scales at a 
lower sampling rate. At the same time, we can also eliminate 
the short time scale components without too much of a 
reduction of the signal quality. The result of this re- 
sampling will enable us to store and transmit a given signal 
with reduced amount of data. 

Data storage is a concern as the HHT is designed to 
examine the details of the signals. For storage, however, 
there are steps HHT can help. These steps are described as 
follows : 

From the designed signals, decompose the data into IMFs . 
Test the recombined IMF without the highest few components to 
determined the degradation of the sound quality. Within the 
tolerable range, eliminate the maximum number of high 
frequency components. Re-sample the data by decimating (i. 
e., with large sampling steps). This will result in a smaller 
number of data, but also with slightly down graded sound 
quality . 

Machine Operating Condition Monitoring 

Referring to Figure 44 the acoustic signal emitted from 
an operating machine can be compared to a kind of speech made 
by the machine trying to inform us of it condition. A chipped 
gear, a broken bearing, a cracked shaft, a worn grinder, 
cutter, or saw blade will all send its distinct sound and 
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particular vibration signals . Such signals can be analyzed 
much the same way as speech signals discussed above. We will 
use a grinder as an example. 

Figure 60a, b and c show the acoustical data from a 
grinder operating on a hard surface continuously for 95, 12 0, 
and 200 hours designed as TR095, TR120 and TR200 respectively. 
From the signals, it is hard to say whether the grinder 
operated at its optimal condition. Of course, the magnitude 
of the signal representing the loudness of the sound is 
increasing, but that alone is not a criterion for determining 
the condition of the grinder, for there might not be a priori 
record of the grinder at other time periods. Also from the 
record, we can determine the rotating speed of the grinder. 
The rotating speed is also changing. For the cases of TR12 0 
and TR200, the rotating speed is 8,000 rpm; while the rotating 
speed for the case of TR095 is slightly higher at around 
14,000 rpm. None of the above characteristics can be used as 
an indication of the health condition of the grinder. 

To determine the health condition of the machine, one has 
to establish a criterion based on the condition of the machine 
at the time of monitoring. This can be implemented by HHT 
analysis. Using HHT analysis, we can first sift the data and 
obtain their IMF components as shown in Figure 61a, b and c. 
Then, we can also construct the Hilbert spectra from the IMFs 
as shown in Figures 62a, b and c. Here some distinct 
characteristics appear. When the grinder is sharp and 
operating smoothly, there is a clear frequency of the acoustic 
signature as in Figure 62a. The energy concentrates around 
the rotating speed as indicated. As the grinder becomes dull, 
it happens first to only part of the grinder surface. The 
noise generated by the dull part causes signal with greater 
amplitude or energy density. This intermittent acoustic 
signature occurs proportional to the dull surface of the 
grinder. This diffusion of acoustical signal can cover a much 
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broader range, but only intermittently. The intermittence, 
however, is only detectable by the HHT analysis as shown in 
Figure 62b. To the ear, the sound might still be continuous. 
As the wearing progresses, all parts of the grinder becomes 
5 dull. At this stage, the acoustical signal will be a diffused 
signal everywhere as shown in Figure 62c. 

The grinder is but only one example. Other telltale 
signs from the chipped gear, the broken bearing, the cracked 
shaft, the worn cutter, or saw blade can be identified in a 
10 similar way. 

To summarize the steps in machine health monitoring, we 

K= 

CJ can state that this task is similar to speech identification. 
l*\ We can establish the templates database of normal and various 
III abnormal operation conditions from the sound and/ or the 
g?| vibration signal of the machine. To diagnosis a machine, just 
record the operational sound and/or the vibration signal of 

s 

□ the machine at any time to make comparison with the templates 

in the database to identify the problems. Knowing problems 
U.I with existing templates, the problem can be identified. For 
r2P safety sake, whenever the machine sound and/or its vibration 
signals are not within the tolerable range of the normal 
condition template, warning signal should be given. 

The dependence on the existence of scale for mode 
25 definition has a limitation in that the decomposition method 
cannot separate signals when their frequencies are too close. 
In this case, there would not be any characteristic scale: 
therefore, physically they are identical. 



3 0 Particular Advantages of The Invention 

The strength of the EMD method should be reiterated. EMD 
is built on the idea of identifying the various scales in the 
data which are quantities of great physical significance. 
Therefore, in the local extrema and curvature extrema Sifting 
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Processes, orthogonality is not a consideration, but scales 
are. Since orthogonal decomposition is a characteristic for 
linear systems, violating this restriction is not a 
shortcoming but a breakthrough. Therefore, the decomposed 
5 IMF's may or may not be orthogonal. As such, this method can 
be applied to nonlinear data. Though the IMF's in most cases 
are practically orthogonal, it is a coincidence rather than a 
requirement of the EMD . 

Another advantage of the method is the effective use of 
10 all the data representing the physical phenomenon. In the 

l s l Sifting Processes, the longest scale is defined by the full 

~~~ 

y length of the data. As a result, EMD can define many long 

ill 

SJ period oscillations. As is well known, the Hilbert transform 

without sifting tends to identify the highest frequency 

tfi (Boashash, 1992, "Estimating- and Interpreting- the 

~ 3 
~- = 

_ 5 Instantaneous Frequency of a Signal, Part I: Fundamentals", 
U Proc . IEEE, 80, 520-538.), the extraction of the long period 

5 U 

components is indeed a new feature of the EMD. 

«; Finally, though the EMD method will give IMF components, 

fib the individual component does not guarantee well-defined 
physical meaning. This is true for all decompositions, 
especially for the methods with a priori basis. In most 
cases, however, the IMF's do carry physical significance. 
Great caution should be exercised in making such attempts. 

25 The rule for interpreting the physical significance of the 
IMF's is that the scales should be clearly separated. 
Together with the Hilbert spectrum, the totality of the 
presentation should give a much more detailed representation 
of the physical processes than conventional methods. 

30 The invention being thus described, it will be obvious 

that the same may be varied in many ways. Such variations are 
not to be regarded as a departure from the spirit and scope of 
the invention, and all such modifications as would be obvious 
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to one skilled in the art are intended to be included within 
the scope of the following claims. 
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